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ABSTRACT 

This report describes a Monte Carlo simulation of transition flow round 
a sphere. Conditions for the simulation correspond to neutral monatomic 
molecules at two altitudes (70 and 75 km) in the D region of the ionosphere. 
Results are presented in the form of density contours, velocity vector plots 
and densitv, velocity and temperature profiles for the two altitudes. Contours 
and density profiles are related to independent Monte Carlo and experimental 
studies, and drag coefficients are calculated and compared with available ex- 
perimental data. The small computer used is a PDP-15 with 16 K of core, and 
a typical run for 75 km requires five iterations, each taking five hours. The 
results are recorded on DECTAPE to be printed when required, and the program 
provides error estimates for any flow-field parameter. 
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1. INTRODUCTION 


1 . 1 Motivation for the Study 

1.1.1 The lower ionosphere 

The term ionosphere is given to the region of weakly ionized gas which 
surrounds the earth from an altitude of 50 km to about 1000 km. 

Ground-based radio techniques have shown that the ionization produced 
by solar and cosmic radiation form a layer structure, which provides a 
convenient method of classifying the various altitude regions. 

Whitten and Poppoff (1965) have designated the region between 50 and 
150 km as the lower ionosphere, and this is subdivided into E and D regions, 
according to the layer structure of the corresponding ionization. 

The E region is a well-defined layer of ionization formed during normal 
daytime conditions in the altitude region between 90 and 160 km. Above this 
is the F region. The D region extends down from 90 km. 

1.1.2 D-region measurements 

Rocket-borne electrostatic probes provide one of the fundamental techniques 
for measuring the properties of the ionosphere. The probe is a small metallic 
electrode carried through the plasma by a sounding rocket. A DC power supply 
in the rocket biases the probe at various voltages positive or negative with 
respect to the plasma and the current collected by the probe provides inform- 
ation about the conditions in the plasma, such as concentrations and energy 
distributions of the charged particles. 

Probe measurements have the unique advantage of being localized, rather 
than averaged over a large volume of plasma, so that the accurate interpre- 
tation of probe data is very desirable. 
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A fundamental aspect of probe flow is the formation of a sheath region 
near the probe, where charge neutrality is not satisfied, as it is in the 
undisturbed plasma. 

Classical Langmuir probe theory is applicable only where the mean free 

path for charged-particle-neutral collision is greater than the sheath thickness, 

so that on the average, particles suffer no collisions after entering the sheath. 

This condition is satisfied only above about 90 km. 

The theory of Gerdien condenser- type instruments treats the motions of 

the charged particles as mobility- control led, which assumes a mean free path 

much shorter than the distance over which the electric field changes 

appreciably. This neglect of space charge effects is not justified for number 

3-3 

densities higher than 10 cm , so that a more general probe theory is required 

13 -3 

for D-region applications, where number densities are greater than 10 cm 

Cicerone and Bowhill (1967) have developed an analytical theory to cover 
a stationary probe immersed in a relatively high pressure, weakly ionized gas. 
1.1.3 Monte Carlo study 

A complete theory of D-region probes must allow for the possible presence 
of two or more negatively charged species and motion of the probe relative to 
the plasma, a problem which at present seems intractable using analytical 
methods. 

A Monte Carlo simulation, made possible by the development of high- 
speed digital computers should be able to incorporate the above features, 
provided that the required collision cross sections are known. The present 
work is concerned with a simpler problem, that of flow about a sphere travel- 
ing. at Mach 2.7 through a stationary monatomic neutral gas with temperature 
and density corresponding to the neutral constituents of the ionosphere at 
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two altitudes (70 and 75 km) in the D region (see Table 1.1). The distribution 
obtained by the present study will be used as a neutral background gas for 
further studies of ion-probe interactions, using the Monte Carlo technique. 

1 . 2 Analytical Treatments of the Shock Problem 

For the purposes of this report a shock wave having supersonic flow 
behind the shock will be termed a "weak shock", while a shock wave having 
subsonic flow behind the shock will be termed a "strong shock". 

1.2.1 The Boltzmann equation 

The state of a gas or gas mixture at a particular instant is completely 
specified for the purposes of analytical kinetic theory if the distribution 
function, f(v, r, t) for the molecular velocities v and position r at time 
t is known throughout the gas . 

Observable properties of the gas may be obtained by suitable averages 
over the distribution. 

The distribution function, f(v, _r, t) is defined such that 

dN = f (v, r, t) dv dr 


is the number of molecules that have velocities between v and v + dv and 
position between r_ and r_ + dr at time t. 

The most general and fundamental description of the time and space rates 
of change of f, due to collisions within the gas is given by the Boltzmann 
equation (Boltzmann, 1964) 


3f 

at 


+ 





coll 
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TABLE 1.1 Free stream data (U.S. Standard Atmosphere, 1962). 
Parameter Height (km) 


_3 

Number density, N (m ) 


70 


1.64 x 10 


21 


75 


.79 x 10 


21 


Temperature, T m ( 0 K) 

215 

195 

Effective velocity, U ro (ms _1 ) 

880 

835 

Molecular mass, m(kg) 

4.81 x 10" 26 

4.81 x 10" 26 

Mean free path, A^Cm) 

1.03 x 10" 3 

2.14 x 10' 3 
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where a(£, t) is the acceleration of a molecule produced by any external 
field when the molecule is at point r at time t, and the right-hand side is 
the Boltzmann interaction term expressing the net rate of change of the 
distribution function at a fixed point due to molecular interactions. 

The form of the interaction term depends on the type of force between 
two interacting molecules. For example for a central force which is a function 
of the distance between two molecules 1 and 2, respectively; 


tdf i /dt) con " / (f i f 2 ' f i f 2 J gb db d£ d — 2 


where f^ = f(vj, r, t) 

f 2 = L> t ) 

f l = L> 

f 2 = f(v.., r, t) 

and (v*, vp are t ^ ie velocities after an interaction of molecules 1 and 2 

which had velocities (v^ v_ 2 ) before interaction, b is the impact parameter 
and e the azimuthal angle of the orbital plane of molecule 2 with respect to 
1 . 

Boltzmann's H-theorem (Boltzmann, 1964) shows that molecular encounters 
will tend to bring about a Maxwellian distribution of velocities if the gas is 
left to itself. 

Thus it may be shown that the Boltzmann equation is satisfied by a 
Maxwellian type distribution function 


f(v, r) = A e 


-C|v| 2 /vh 


where A is a constant and v is the most probable thermal speed of molecules, 

m 


given by (2kT/m) 


1/2 
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Boltzmann's H is defined by 

H(t) = /“„ f(v, t) In f(v, t) dv . 

A minimum of this function is a necessary condition for equilibrium and 
it may be shown that the above Maxwellian distribution for f satisfies this 
condition. 

Efforts to obtain exact solutions to the Boltzmann equation for the 
more general case of non -equilibrium have proved much less successful, both 
because of the non-linear nature of the Boltzmann equation and because of 
the intractable form of the collision integral. 

The most promising analytical approach seems to be the assumption of 
an explicit form for f in each particular case (Mott-Smith, 1951), see 
Section 1.2.3. 

1.2.2 Analytical solutions 

The equations of fluid flow may be obtained by solving the Boltzmann 
equation for the space and velocity distribution of the molecules by the 
Enskog- Chapman method (Chapman and Cowling, 1939). 

Boltzmann's equation is expressed in the form 5(f) = 0 and a series 
solution 

f-i f tr) 

r=o 
( rl 

is assumed, where f is the rth order term in the expansion of f. 

This leads to a solution for the distribution function in terms of 
the parameter X/Ax where X is the mean free path and Ax the distance in 
which f changes by an appreciable fraction of itself. 
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The zero-order terms give the equations of flow of an inviscid fluid. 

The first-order terms give the Navier-Stokes equations, and the second 
order terms, the Burnett equations. The simplest inviscid flow problem treats 
a one dimensional flow of gas produced by the motion of a piston in the axial 
direction, within a cylindrical tube (Becker, 1968). 

In this case the equation describing the motion reduces to 


P 



dus 

u 33? 


+ c 


2 o 

s dx 


( 1 . 1 ) 


where p is the density, u the velocity at a point, and C g is the local speed 
of sound. 

The method of characteristics introduces two families of curves t^. 


dx 

dt 


u + 


C 

s 


dx 

dt 


u - 


C 

s 


which are characterized by parameters X and y respectively, then Equation (1.1) 
and the continuity equation may be expressed in a (X, y) coordinate system 
as 


where 


w + u - const on t ^ 
to - u — const on t ^ 

oj(p) = J — dp 
J P 0 P 


( 1 - 2 ) 


The two relations (1.2) then form the starting point for computing p(x, t) 
and u(x, t) in a specific problem, i.e., with given initial values. 
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Solutions of the Navicr-Stoke equation were obtained by Becker (1923), 
assuming constant coefficients of viscosity and thermal conductivity. These 
solutions were improved by Thomas (1944). 

For strong shocks the thicknesses calculated by Becker and Thomas are 
very different, but are at most of the o/’er of a few mean free paths. This 
throws doubt on the valiuity of the Navier-Stokes equations , since they are 
valid only if f changes only by a small fractional amount in a mean free path. 

1.2.3 The Mott -Smith method 

In the Enskog- Chapman theory, f is represented by a skewed Maxwellian 
form, having only one strong maximum. Mott-Smith (1951) has suggested that 
a more profitable assumption might be a bi -modal form. Here the distribution 
is assumed to be the sum of two Maxwellian terms (representing subsonic and 
supersonic streams) with different temperatures and mean velocities but with 
unassigned space densities. The densities re obtained from the solution of 
a transport equation for u n , where n is an integer and u is the component of 
molecular velocity in the stream direction. 

Since Mott-Smith's method does not take account of interactions between 
particles of the same stream, his theory is most suitable for strong shocks 
where these are less significant. 

The two- fluid model has been improved by Ziering, et al. (1961), and 
their results are in good agreement with experiments performed by Sherman 
and Talbot for both large and small Mach numbers. (Mach number, M = local 
flow speed/ local sound speed) . 

However, the Mott-Smith models predict the wrong value of the Prandtl 

number, Pr * C h/k (where C = specific heat at constant pressure, n = coef- 
P P 

ficient of viscosity, k = thermal conductivity) near the downstream boundary, 



so that a completely satisfactory analytical solution for arbitrary shock 
strength is still lacking. 

1 . 3 Monte Carlo Evaluation of the Boltzmann Collision Integral 
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Nordsieck and Hicks (1966) have devised a Monte Carlo method for the 
evaluation of the Boltzmann collision integral. The method has a major 
advantage over the above analytical approaches in that it can easily be 
modified to use any molecular force model as long as the differential cross 
sections are known. It may also be used to test any velocity distribution 
function proposed as an approximate solution of the Boltzmann equation for 
a shock wave or other flow conditions, or to check directly the various 
elaborate analytical calculations involved in moment methods (Martikan, 1966). 

The procedure replaces the collision integral by an integral over a 
finite region of velocity space, taken so as to include most molecules. The 
average of the integrand over all values of the line of centers vector is 
then approximated by the average of a large and fair sample of particular 
values of the integrand, selected by Monte Carlo trials. 

The method has been applied to both the pseudo shock (a translational 
relaxation of molecules) and the shock structure. 

Extensive error analyses performed by Hicks (1968) have estimated that 
for a Mach number of 2.5 the random errors in the velocity distribution function 
and the collision integral amount to 2% or less, and random errors in moments 
of these functions range from 0.03 to 2.7%. The complete program required 
8,000 words of storage on a CDC 1604 computer. 

A major difficulty with the practical application of approaches based on 
the Boltzmann equation is the inclusion of realistic and complicated boundary 
conditions, especially surface interactions with bodies placed in ths flow. 
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In an attempt to solve problems of this type, several Monte Carlo algorithms 
have been developed which treat the dynamics of gas molecules more directly. 

1.4 Monte Carlo Simulations of Rarefied Gas Flow 

Bird (1965) has developed a Monte Carlo technique which was applied to 
the problem of a gas initially in equilibrium between two infinite, plane, 
parallel and specularly reflecting walls. One wall then impulsively acquires 
a uniform velocity towards the other, and the numerical experiment studies the 
shock wave so formed. 

The first step in the procedure is to select a molecule at random, thei; 
sample the number density in the vicinity of the molecule, which is retained 
or rejected such that the probability of retention is proportional to the local 
density. A second molecule is chosen at random, subject to the condition that 
its position be within half a local mean free pa '' of the first molecule, on 
the right side for a collision to occur. The relative velocity is determined, 
and the pair is accepted or rejected so that the probability of retention is 
proportional to the relative velocity. When a pair is retained, a line of 
impact is chosen at random and new sets of velocity components for the 
molecules are computed. Each time such a collision takes place, a new pair 
is selected and the time is advanced by 

At .1- 

v 

where N is the number of molecules used in the simulation, and v is the 
o 

collision frequency from kinetic theory. After a suitable time, t^, all the 
molecules and the wall are moved through a distance appropriate to t and their 
current velocities. At much larger time intervals the velocity and density 
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profiles between the walls are sampled. After a run of thirt> minutes on the 
Atlas digital computer at the University of Manchester, England, Bird obtained 
shock profiles with a standard deviation of two to four percent. The shock 
was found to travel at the speed predicted by the Rankine-Mugoniot equations 
and resulting profiles showed good agreement with the Mote-Smith results, at 
a shock Mach number of 1.5. 

Vogenitz, et al . (1968) have applied Bird's method to transition flow 
about cylinders, spheres, wedges and cones, the results were compared with 
wind tunnel tests using a free molecular recovery temperature probe, and with 
electron beam density measurements. Agreement in both cases was good. The 
computation times ranged from 5 minutes using 32,000 words of storage to 20 
minutes using 80,000 words, on a CDC 6600 machine. 

Bird's technique has been used to provide pictorial simulations of 
transition flow (Bird, i969j , which should prove to be a valuable visualization 
technique for this difficult area. 

The present program of study is intended to show that comparable results 
may be obtained using much smaller compute^ capacity. The machine used is a 
Digital Equipment Corporation PDP15, having 16,000 words, each of 18 bits. 

The machine is intended primarily for on-line data reduction of ionospheric 
observations of the partial reflection type and the Monte Carlo simulation 
can be run whenever the machine is not required for this work, which takes 
place only during daylight hours. Monte Carlo studies might well prove to 
be valuable users of the small on-line computers having relatively low duty 
cycles, which are to be found in many installations about the country. 
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2 . THE PROGRAM 


2.1 Monte Carlo Choice 

The core of any Monte Carlo method is the Monte Carlo choice. TTie process 
to be modeled is reduced to a series of decision points, which are branch 
points in the process, where any one of a number of future courses is possible. 
Probabilities must be assigned to the possible events, according to some theory 
of the microscopic kinetics of the system under consideration. The probabilities 
are mapped onto the interval 0 to 1 as shewn in Figure 2.1, so that the fraction 
of the unit interval allotted to each event is equal to its probability. 

A random number is now chosen from a distribution uniform between 0 and 
1. The event is determined by placing the random number on the probaoility 
interval; if it falls in the interval allotted to event A, then event A is 
said to have occurred, if in the interval of S, then B occurred and so on. 

Since the random number may fall with equal probability at any point on the 
unit interval, the probability of its falling in A is exactly equal to that 
fraction of the unit interval allotted to event A, namely the computed prob- 
ability of event A. In this way events in the real process are modeled on 
the computer, provided only that a random number generator of fairly uniform 
distribution is available, and that realistic probabilities can be assigned 
to each possible event. 

In the present simulation it is also necessary to draw particular events 
from a continuous distribution, such as initial position and velocity of a 
molecule or velocity after a collision. In this case a modification of the 
above procedure is used. The following Equation (2.1) from probability theory 
relates a random variable x, with probability density function f(x), to a 
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Figure 2.1 Monte Carlo choice. 
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random variable uniformly distributed on the interval 0 to 1 (Shreider, 1966) 


R. 

l 


X f(0 dc 


x 


1 


( 2 . 1 ) 


where Z is a dummy variable 

Xj is the lower limit of the particular 
x in question. 

Using this relation and the known distribution for the parameter to be 
drawn, the Equation (2.1) is integrated and rearranged to give x as an explicit 
function of . A random, or as in the present program, a pseudorandom number 

is then drawn, and used to e/aluate the corresponding x. In this way if the 

R^ values are distributed between 0 and 1, with distribution approximately 
uniform for a large sample, the resulting values of x will have a distribution 
close to f (x) . 

2.2 The Method 

2.2.1 Characteristics of the gas 

In the present simulation of rarefied gas flow, the real events are 
collisions between the molecules in the region about the probe. 

Using an approach analagous to that of Mott-Smith, the gas near the probe 

is considered to consist of three distinct classes of molecules, each character- 
ized by a distribution of density, mean velocity and temperature. 

Class 1 molecules have not yet encountered either the probe, or a 
member of any other class . 

Class 2 molecules have been reflected specularly from the probe, 
and have only encountered other class 2 molecules since reflection. 
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Class 3 molecules have encountered at least one molecule of another 
class, either before or after striking the probe. 

The method involves an iterative procedure which gives successively closer 
approximations to the actual distributions for the three classes. 

2.2.2 Division into cells 

The continuous distributions of velocity and temperature are approximated 
by dividing the space of interest into discrete cells, each having one value 
of every parameter for the three classes of molecules. 

Since the probe has axial symmetry, a cylindrical coordinate system is 
used (Figure 2.2) having ten cells in the radial and 30 in the axial directions. 

The resulting system of 300 cells, are coaxial cylindrical shells (Figure 
2.3) and it is assumed that the flow has no mean variation in the azimuthal 
direction. 

Since class 2 and 3 molecules are found to be confined to the region of 
the probe, class 2 velocities are stored only for the ten axial cells nearest 
the probe, and class 3 velocities and temperatures only for the 20 axial cells 
nearest the probe. 

Beyond these regions class 2 mean velocities and class 3 mean velocities 
and temperatures are set to the last computed value for that radial shell. 

2.2.3 Computational procedure 

This section gives a brief outline of the whole computational procedure, 
which will be dealt with in more detail in sections 2.3 and 2.4. 

A first approximation to the various distributions is taken as a back- 
ground gas, into which are introduced test molecules of class 1, starting from 
the entrance plane. 
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A test molecule is introduced having position and velocity coordinates 
chosen at random from distributions calculated for the flux across the input 
plane. 

The molecule is examined after a fixed interval of time 6t, when the 
Monte Carlo choice probabilities are calculated, and a random selection is 
made to determine whether a collision has occurred in the last time interval 
fit, and if so the class of the collision partner. 

If a collision has occurred, new velocity components are selected, accord- 
ing to the laws of classical gas dynamics. 

The velocities and squared velocities of the molecule are recorded for 
that position, and a count of particles for the test class is incremented 
by one. The test molecule proceeds by straight line path segments towards the 
probe, and at each point the position is tested to determine whether the test 
molecule has struck the probe, or reached the boundaries assigned to the region 
of the study. 

If a class 1 molecule reaches the radial boundary it is reflected specu- 
larly, so as to model the introduction of new molecules from an infinite real 
gas. If a test molecule reaches either axial boundary without striking the 
probe, a new test molecule is introduced at the entrance plane and the whole 
process is repeated. 

If a class 1 molecule strikes the probe, it is reflected specularly and 
becomes a class 2 molecule. The mean temperature is assumed to remain unchanged 
after reflection. 

After a suitably large number of molecules have been followed in this 
way, the accumulated parameters are used to calculate a new approximation to 
the background gas distributions. These values replace the old background 
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gas for the next iteration, when the introduction of test particles begins 
again. 

After a number of iterations, the distributions converge to a stable form 
with steady statistical fluctuations, which is taken to be the final estimate 
of the model flow. 

2 . 3 Description of the Program 

The flowchart. Figure 2.4 shows the logic of the Monte Carlo program. 
Subroutine names appear in wide-spaced letters in Figure 2.4 and capitalized 
in the text. A complete listing of all programs appears in the Appendix. 

2.3.1 Cell system and storage of parameters 

As explained in section 2.2.2, the cell system consits of ten cylindrical 

shells, each divided axially into 30 equal rings, generating 300 cells so that 

300 values of each parameter must be stored (with the exceptions noted in 2.2.2). 

In the present study, core storage is limited to 16,000 18-bit words, so 

that economy in storage is an important criterion. 

Each parameter is stored as a two-subscript array, with the first subscript 

representing axial, the second radial cell numbers. The two subscripts for 

a given position are found by taking each coordinate, dividing by the corresponding 

cell dimension and rounding to the next higher integer. 

The PDP-15 requires two 18-bit words to store a real number, as shown in 

Figure 2.5. This provides an accuracy of six decimal digits. An integer con- 

17 

stant is stored in one word, giving a maximum magnitude of 131071 or 2 -1. 

It is felt that integer storage can provide adequate accuracy for all background 

gas parameters, while real-variable storage is used for all accumulated velocity 

parameters to allow for the large magnitudes involved. Number densities are 

-17 

normalized before storage by multiplying by a constant factor (F.. ) of 10 



Figure 2.4 Main program flow chart 
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17 

The reciprocal CF j)» 10 is used to restore the true density before use 
in the program. 

2.3.2 Read parameters, initiali z e subroutine RANDU, and arrays, 
subroutine SETUP 

At the beginning of any run, and only then, the parameters defining the 
conditions of the simulation are read from the teletype keyboard. These 
define for example, cell size, ambient free stream conditions and probe posi- 
tion (see Section 2.6). 

A prime integer is supplied to initalize the random number generator 
(Section 2.4.1) and the subrountine SETUP sets fixed parameters (such as it) 
and initial values of all arrays, to free stream conditions for class 1 back- 
ground parameters, and zero for all accumulated parameters. 

2.3.3 Introduction of a new molecule, subroutines NEWPOS and NEWV 

The introduction of a new molecule involves two subrountines NEWPOS 

and NEWV which are explained in more detail in sections 2.4.2 and 2.4.3. 

NEWPOS generates the cylindrical coordinates of the initial position 
(0, r, 0) on the entrance plane. The r and 6 coordinates are chosen from dis- 
tributions uniform with area over the entrance plane. 

The cell number subscript is assigned as in section 2.3.1 and NEWV gen- 
erates cylindrical coordinates of the new initial velocity (V z , V^, 0 ) for 
a flux of molecules with mean axial speed U^. (Figure 2.6) 

2.3.4 Tracing a molecule, subroutine INCPJS 

Basic to the operation of the program is the small time increment 6t 
at which all parameters of the test molecule are recalculated; this time 
increment will be described hereafter as a "mo”. 

The time is incremented by the chosen 6t, and subrountine INCPOS returns 
the next position using the molecular velocity and previous position coordinates. 
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The molecule is assumed to traverse a stra' ght line path between collisions. 

2.3.5 Tests for cell and system boundaries, subroutines NEWCEL, HITSPH. 

RBOUND 

The next operation is to make a series of tests on the position coordi- 
nates of thw molecule, to determine whether it has moved into an adjacent 
cell, and if so whether that cell lies outside the chosen system boundary, 
which is drawn so as to include the surface of the sphere (sec Figure 2.2). 

If a class 1 molecule has crossed the radial system boundary, the sub- 
routine RBOUND simulates the introduction of a new class 1 molecule from 
the gas beyond the system, by means of a specular reflection. Class 2 or 

7 ? . ■ . 1 j*_i u «... r J « 1 ~ 4- ^ 4 r r ^ 
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exists from either plane axial boundary (the entrance or exit plane), that 
molecule is lost and a new molecule is chosen at the entrance plane. 

If the molecule is found to have struck the probe it is reflected specu- 
larly. Subroutine HITSPH calculates new position and velocity components 
using the probe geometry (see Section 2.4.7). 

2.3.6 Collision probability, subroutine VREL 

If the molecule has not crossed a cell boundary in the last mo, values 
of collision frequency with the three background classes (v^, v^, v^) and total 
collision probability, P(bang), remain the same as for the previous mo. Other- 
wise new collision frequencies are generated using the formula, 


where v. is the collision freauencv with background class i 
i ‘ 

is the number density of background class i 



V . is the mean effective collision speed for the test and ith 
reli 

background class, given by 


GO 00 00 



00 OC 00 


where v ^ is a velocity of the ith background class, having Cartesian 

components v biy v biz an< ^ a ^ axwe ^-*- an distribution 

function P(v , . ) . 
v is the test molecule velocity. 

This integral is evaluated in Section 2.4.4, and the function VREL returns 
a value of V using a 15 step approximation to the integral. 

2.3.7 Collision leading to class three test subroutines MONTE, BANG, 
PTNR3 , VELS 

A Monte Carlo choice (see Section 2.ij now determines whether a collision 
occurred in the last mo. RANDU returns a random number which is compared with 
the collision probability P(bang). If it is less than P(bang), a collision 
is said to have occurred and the subroutine MONTE makes a further Monte Carlo 
choice tc determine c^ the class of the collision partner. Collision proba- 
bilities for each class are computed from the corresponding frequencies (see 
Section 2.3.6), and mapped onto the unit interval. A random number then deter 
mines the choice of c^. 

Subroutine BANG, determines the combination of test class c and partner 
class Op, so that if the two are different or c^ = 3, the correct combination 
of parameters is passed to suboutine PTNR3. This subroutine determines the 
test molecule velocity after collision, using a hard sphere model and incorpor 
ating the mean persistence of velocities after collision. Subroutine VELS 
provides a random thermal velocity, drawn from a Maxwellian distribution with 
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spherical symmetry as in Sections 2.4.5 and 2.4.6. The test molecule now 
moves into class 3. If the test and partner molecules were found to be of the 
same class, the collision is ignored, unless both were of class 3, in which 
case the collision is computed as described above. This ensures that the pro- 
cess remains collision dominated even after class 3 molecules predominate near 
the probe. 

2.3.8 Bookkeeping and iterative procedure, subroutines BKKEEP, NEWBG 

At this point, subroutine BKKEEP increments the accumulated parameters 

for the present class of test molecule, i, at its present position as determined 

by the array subscripts. A count of molecules Nh is incremented by 1, and for class 

2 or 3 molecules the current velocity is added to accumulated velocity sums 

IV . , IV . for the axial and radial directions. For class 3 molecules only, the 
zi ri 

2 2 

squared velocities are cdded to accumulated squared velocity sums I(V_ 7 ) , I(V r7 ) 
which will be used to calculate mean class 3 temperatures. Again the velocity ac- 
cumulations occur only in the region of the probe, as explained in Section 2.2.2. 


The procedure then begins again for the next mo. 

Each test molecule is traced and recorded until it leaves the system 

through one of the boundaries. A new class 1 test molecule is then introduced 
at the entrance plane, and the tracing process is repeated. 

After a suitable number of molecules have been traced, as determined by 
standard deviation estimates (Section 3.1), a new background gas is computed 
from the accumulated parameters for each class in each cell as follows. The 
new number density NL is given by 


N. = 

l 


M. x F 
l norm 

V cell 


where 


M. is the accumulated number of molecules of class i recorded in 

l 


that cell 
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by 


v is the cell volume 
cell 

F is a constant normalizing factor for the system such that the 
norm 

total number density at the cell farthest from the probe is that 
of the freestream. 

The new mean velocity components for classes 2 and 3, (V are given 


I v • 

r, L Zl 


V . = 
zi 


M. 

l 


( 2 . 2 ) 


n 


L rl 

M. 

l 


(2.3) 


where Tv . , 7v . are the sums of ail axial and radial velocities respectively 
L zi’ L n 

recorded in that cell as explained in Section 2.2.3. 

The new class 3 temperatures (T z3> T r3 ) are £i ven by 






(m 3 -D 




( M 3-D 


(2.4) 


(2.5) 


where m is the mass of molecule 
k is Boltzmann's constant 

9 9 

[(V ) , l(V t) are the sums of squares of all axial and radial 

ZO 1J 

velocities recorded in that cell (see Section 2.2.3) 

V ,, V , are the mean velocities of classes 2 and 3 as above 
z3* r3 

M 3 is the total number of class 3 molecules recorded at that cell. 


Class 1 temperatures and velocities and class 2 temperatures are assumed in- 
variant with position everywhere. Class 2 and 3 velocities and class 3 tem- 
peratures are stored only for the region nearest, the probe as explained in 


Section 2.2.2. 



/ 


28 

After the new background parameters have been calculated, all accumulated 
parameters are reset to zero and the new iteration begins with the introduction 
of the first molecule at the entrance plane. 

2.3.9 Program output, subroutines RESULT, GIVE 

The subroutine RESULT can take various forms depending on the output re- 
quired. The standard form, included in the binary library file (.LIBR5 BIN) 
assembled for this program is called RESGIV. This uses subroutine GIVE (not 
shown on the main flowchart for clarity) to write the accumulated parameters of 
each iteration on DECTAPE in a non-file oriented mode, using the DTF DECTAPE 
handler. Each WRITE command fills 256.^ (256 decimal) words (one block) of 
the tape from 256^ words of core. Unused words within this number are filled 
with blanks. The data to be written have been arranged so as to fill the tape 
in an economical fashion using this mode. Since unc DECTAPE comprises 576.^ 
blocks, up to 36 iterations, each of 16 blocks may be stored in this way. 

Alternatively RESNM or RESRAW may be loaded in place of RESGIV, to output either 
raw data (RESRAW) or data normalized by free stream values (RESNM) on the tele- 
type. Which of the 1900 parameters are printed each time is of course a matter 
of choice. 

Core storage limitations permit the loading of only one form of RESULT 
for any given run. 

2.3.10 Off-line data reduction, subroutines NEWBG5, RESNM7, STDV7, RSTDV7 

When the required number of iterations have been recorded on DECTAPE, 

using RESGIV and GIVE, execution of the main programs may be terminated, and 
one of several possible sets of data reduction program loaded instead. NEWBG5 
and RESNM7 read the accumulated parameters of each iteration from DECTAPE, com- 
pute the corresponding background gas and print the normalized background 
parameters on the teletype. 

i 
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STDV7 and RSTDV7 may be used to process the short runs (see Section 3.1) 
to give estimates of the mean and standard deviation for a full scale run. 

2.3.11 Dumping the program, subroutines DUMP, LSSW, PAWSE 

The PDP-15 single-user Advanced Monitor Systen KM15 V4A allows for the 
possibility of writing the entire core onto DECTAPE at any time during execu- 
tion. The process is initiated by pressing the keys CTRL and Q simultaneously 
on the teletype, and is called a dump. An area of DECTAPE must be reserved for 
dump, and the dump may be recovered and execution restarted exactly as if no 
interruption occurred, provided the program is in a suitable waiting loop such 
as a Fortran PAUSE or the Macro subroutine PAWSE used here, when the dump is made. 

The Macro routine LSSW provides a way of initiating a recoverable dump 
at any time by depressing a specified console data switch. The subroutine 
DUMP is then called and the dump may be initiated using the CTRL and Q keys. 

After recovery of a dump, DUMP also repositions the data tape correctly, so that 
the entire computer and all tape drives are available for other uses between 
dumps. In this way execution of the program may be resumed each night, while 
the computer is otherwise occupied during the day. 

Subroutine LSSW proved very valuable during modification and de-bugging 
of the program, since for example, the position of a particle can be printed 
when required by the user, rather than every time the program reaches a given 
point (see the program comments in Appendix). 

The Macro subroutine PAWSE may be used instead of a Fortran PAUSE state- 
ment, causing a halt in execution until the keys CTRL/P are struck on the 
teletype. The Macro program is considerably smaller than the corresponding 
Fortran routine since it cannot handle numbered PAUSE statements. This 
facility is not required for the present purpose, and PAWSE is used here to 


save core storage. 



2.4 Important Subroutines 
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2.4.1 Pseudorandom number generator, subroutine KANDll 
The pseudorandom number generator used here is of the multiplicative 
congruent type, based on CACM algorithm 294 (Stromc, 1966). 

The program generates the next uniformly distributed pseudorandom number 
on the interval (0, 1) as in the flow chart (Figure 2.7). The procedure uses 
two constants, M and C, chosen to maximize the period and minimize the correlation 
of the sequence generated. For the present program the following equations 
are used for choosing M and C 


M = D 

where D = number base of the machine 

k = entier ((2n+l)/3) i.e., the integer part of ((2n+l)/5) 
where n is the maximum number of significant digits for a real variable stored 
in the machine. 

Then . 

r n n " k 

C = D - q 


where 


q = 3 for D n " k > 100 


( 2 . 6 ) 


Fortran double precision arithmetic (8 significant decimal digits) 
yields; D = 10, n = 8 

k = entier (2x8+ l)/3 = 5 
n-k = 3 

i.e. M = 10 5 

C = 10 3 - 


3 = 997 
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DU should be a positive integer less than and rclathely prime to M, here 
I 'U was chosen to be 8287. 

Difficulty was experienced in trying to use si v.lo precision arithmetic 
only, for a pseudorandom number generator or. the POP- 15. 

While six significant decimal digits are claimed bv DEC for single 

precision real arithmetic, the sixth decimal digit did not prove reliable 

enough for this use, so that after a series of calls of the program, incorrect. 

digits began to appear in the lowest significant decimal places. These were 

increased in significance by successive multiplications, giving the sequence 

of pseudorandom numbers undesirable properites. The effect was a non-uniform 

distribution of numbers selected at large, random intervals from the sequence, 

which appeared as an incorrect distribution of samples with radius at the 

entrance plane. Assuming only five significant decimal digits leads to a 

n — k 

violation of the requirement (2.6) as D = 100. 

The double precision version used in the program assumes only eight 
significant decimal digits of the nine claimed by DEC. 

Figure 2.8 shows the results obtained with the double precision version 
of RANDU. Numbers introduced in each radial cell N are plotted against 2ir-l 
where ir is the cell number. A perfectly uniform distribution with area 
would result in the straight line shown, with the standard deviation limits 
predicted by a Poisson distribution. The actual points are seen to be in good 
agreement with this hypothesis. 

2.4.2 Selection of position at input plane, subroutine NEWPOS 

The distribution of particles over the input (z = 0) plane must be uniform 
with area, as in the real flow. 
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New z and Q coordinates may be chosen at once; 


2 = 0 


e , 2, R 0 

where R is a random number chosen from a distrihn i >n uniform on the interval 

u 

0 to 1. 

The distribution function, f , for r is given ly 

f -= Q — 

r 2 
r 

max 

where r max is the radius of the system, so that the probability of selection, 
P(r), within any radius r, is given by 


2r 

— dr' 


a fraction proportional to the area of the disc, radius r. 

The distribution is normalized so that the total probability of selection 

within the cell svstem, P(r ) is 1. 

max' 


nr max 2r' . , . 

P(r ) = — = — dr' = 1 

max' 2 

' o r 

max 


Setting this distribution function in the general relation (2.1) gives 


r 2 


dr' = 


'or r 

max max 


where is a random number uniformly distributed on 0 to 1. 


i . e . r = r >^ / R 

max r 

This result is then used to select values of r, using random numbers R f . 
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2.4.3 S election of velocity at entrance plane, subroutine NEWV 

Consider the velocity components of the input class 1 molecules drawn 

from the flux crossing the imaginary input plane, at a mean axial speed U^. 

The class 1 molecules are assumed to have a Maxwellian velocity distribution 

with temperature T^, so that the number of molecules crossing the input plane 

in the +z direction, which have thermal velocity components in the range 

dv , dv , dv about v , v , v is 
x y z x y z 

7 2 2 2 

N -(v" + v + V )/v 

dq = ■ (v U ) e X y 1 m dv dv dv 

n 3/2 3 z •» x y z 

7T V 

m 

where is the free stream number density 

v^ is the most probable thermal velocity in the free stream given 
by (2kT 1 /m) 

The total flux q of molecules crossing the input plane is found by inte- 
grating this expression for v^ and v^ frcm -°° to +°° and v^ from -U^ to +», 
the latter limitation excludes particles crossing in the -z direction (Fan, 1967) 


giving 


,cc 

J 

8 

V =U ■' 

J 


2 Oo \j =-co V tr-oo 

y x 


dq 


q = 


2r. 


1/2 


m 


X(S 2 ) 


where x( s 2 ) is a function of the flow velocity in free stream 


X(S Z ) = e Z + n 1/2 S z [l+erf(S z )] 
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where S = U /v 
z 00 m 

The distribution function of molecules crossing the surface is given 


by 


= 2 

q it 


1 



( 


+ S ) e 
v z 
m 


2 2 2 

-(v +v +v )/v“ 

x y 2 m dv d V d V 

x y : 


which represents the number of molecules with thermal velocities in the range 

dv , dv , dv about v , v . v as a fraction of the total number crossing 
x y z x y z 

ihe input plane in the + z direction. 

This is best re-written in the cylindrical coordinates of the present 

study; V , V , 6 
1 z r v 


dq 2 1 3 , z c ■. 

— 1 77; — r- v ( — + S )e 

q tt x(Sj m v 


m 


2 2 2 

-(V ‘ +V )/v 

r z "V dv de dv 

r r r z 


(2.9) 


The marginal distribution functions for the three components V , V^, 

6 , are each obtained by integrating over the ranges of the remaining two 
r 


components . 
Thus 


/■“ r2iT ? -V 2 /v 2 

f = V f de dV = =-- V e r m 

r J r v z 2 r 

i, o v 

U m 

z 


- f 

6v J 

-U 


V fdV dv 
r r z 


1_ 

2tt 


and 


f v = l 


2 TT 


V fdV de = v ■ - 

r r v x(S )v V v z 
* • z m m 


V -V 2 /v 2 
, z^ c ^ z m 
( — ►S )e 


0 
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The Equation (2.1) relates the above distributions to a uniform distri- 
bution on the interval 0 to 1. 

If the random variable x has a probability density function f(x), then 
the distribution of the random number R. is uniform on the interval 0 to 1 

l 

where 


R. 

i 


f x 

fU) d£ 


where x^^ is the lower limit of x. 

Applying this to the velocity distribution functions for the input 
test molecules gives 


R 


r 


e 



2 

m 


R 


e 


a 

v 

2tt 


and 


R 


x(s z ) 


-(v/v r 

z m 


+ir^^S [1-erf (V /v ) ] 
z l v z m J 


( 2 . 10 ) 


where R , R^ , R^ are independent uniformly distributed random numbers on 
the interval 0 to 1. 

In order to select velocities, the above three equations must be re-arranged 
so that V r> 0 v> and appear explicitly, in terms of R^ R Q R^ respectively. 

In the first two cases this is simple, giving 


V 

r 


v 


m 


/RT/fT) 
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and 

0 = 2tt R 

v e 


However the z component of thermal velocity cannot be expressed explicitly 
in terms of R^. Instead the Equation (2.10) is fitted by a polynomial, using 
a standard computer routine. 


For a given value of S z , 50 pairs of corresponding values of (V^/v^) 
and are computed. It was found that good results could be obtained by 
a least squares fit to a polynomial of seventh degree in terms of n where 


n = |ln(l/R z )| 1/4 . 

So that the corresponding 50 values of n ana V /v are fed into the 

z m 

subroutine, yielding the coefficients a Q through a^ for the equation 

,, , 2 3 4 5 6 7 

V /v = a^ + a n+a 0 n +a n +a.n +a c n +a,n +a_n 
zmol23456 7 

The coefficients are tabulated in Table 2.1 for 14 values of S in the 

z 

range of interest (Fan, 1967). 

Once the coefficients are known for a given S , values of R can be fed 

b z z 

into the polynomial yielding any number of velocity selections. 

2.4.4 Collision probabilities, subroutine VREL 

The following is a more detailed explanation of the theory used in choosing 
collision probabilities. The argument follow Jeans (1954). 

Consider a test particle of speed c colliding with a molecule chosen 
from one class of the background gas, of number density n, having speed c* . 



TABLE 2.1 Polynomial coefficients for subroutine NEWV 


39 



rH 

US 

sO 

to 

o 


LO 

oo 

as 

c- 

sO 

as 

LO 

to 


LO 

fO 

o 

C- 

as 

sO 

00 

as 

as 

as 

CM 

rH 

r^> 

r-. 



to 

to 

to 

o 

rH 

o 

to 

rH 

rsi 

C- 

00 

sO 

LO 

r*". 

o 

00 

to 

Tf 

00 

oo 

os 

to 

o 

to 

CM 

rH 


to 


o 

o 

to 

00 

oo 

sO 

sO 

c- 

© 

rH 

o 

os 

rH 

00 


o 

1 

o 

o 

© 

o 

o 

1 

o 

1 

o 

o 

o 

o 

1 

o 

1 

rH 

o 

1 



to 

to 

o 

rH 

CM 

o 


o 



00 

to 

f- 


r—i 

LO 

so 

00 

sO 

o 


i-H 

00 

00 

00 

00 

sO 

CM 


US 

sO 

to 

as 

LO 

oo 


sO 

sO 

CM 


as 

CM 

00 

sO 

r>. 


so 

o 

o 

r- 


sO 

rH 

oo 

*•3* 

as 

r- 

CM 

aJ 

o 


to 

as 

CM 

CM 

to 

LO 

© 

s O 

sO 

CM 

LO 

os 


o 

i 

o 

1 

rH 

1 

to 


to 

to 

to 

o 

1 

o 

O 

'O* 

LO 

to 


00 

00 

VO 

to 

CM 

LO 

to 

o 

•sj- 

rH 

CM 

LO 

to 

to 


LO 

SO 

Os 

i-H 

rH 

sO 

CM 

LO 

to 

to 

sO 

as 


00 

LO 

o 


00 

to 

rH 

’M' 

sO 

rH 

so 

© 

sO 

LO 

as 

as 

rt 

CM 

00 

O 

o 

i-H 

o 

rH 

rH 

LO 

LO 


00 

o 

CM 


CM 


OS 

CM 

o 

sO 

00 

CM 

o 

00 

sO 

00 

sO 

to 


o 

o 

O 

LO 

1 

vO 

i 

Tf 

1 


LO 

1 

o 

o 

CM 

LO 

00 

LO 

1 


00 

as 

OS 

sO 


o 

os 

to 

as 

o 

as 

as 


rH 


to 


CM 

to 

to 

to 

o 

to 

rH 

oo 

OS 

LO 


i-H 


00 

LO 


sO 

rH 

to 


to 

•** 

LO 

to 

00 


LO 

rj- 

rH 

CM 

rH 

CM 

s0 

c- 

to 

so 

as 

c- 

SO 

c- 

sO 

rH 

rt 

i-H 

to 

CM 

o 

rH 

LO 

CM 

as 

o 

rH 

00 


as 

o 


o 

1 

o 

to 

CM 

rH 

i-H 

1 

l-H 

1 

o 

1 

o 

1 

rH 

to 

rH 

1 

o 

CM 

1 


s£> 

CM 

OS 


rH 

LO 

rH 

rH 

to 

Ht 

sO 

o 

as 

to 


1 — < 

CM 

i-H 

rH 

sO 

SO 



CM 

CM 


as 

LO 



sO 

LO 

rH 

sO 

r-. 

sO 


00 

as 

SO 

LO 

os 

00 

CM 

tO 

sD 

00 


LO 

o 

CM 

CM 

CM 

r- 

as 

CM 

r- 

r- 

00 

cd 


CM 

oo 


r- 

oo 

LO 

LO 

o 


os 

to 

rH 

CM 


o 

Cs| 

LO 

rH 

rH 

o 

o 

o 

o 

to 


CM 

rH 

CM 


1 

1 

1 

rH 

rH 

rH 

rH 

rH 


1 

1 

i— » 

i-H 

rH 


LO 


o 

oo 

LO 

to 

rH 

t- 

as 


CM 

sO 

as 



o 

to 

c- 

s O 

to 

SO 

os 

sO 

rH 

to 

oo 

00 

LO 

CM 


CM 

CM 

CM 

LO 

CM 

CM 

o 

sO 

rH 

os 

CM 

vO 

to 

o 


o 

to 

lO 

LO 

as 

SO 

rH 

as 

C- 

o 

i-H 

as 


as 

CM 

c- 

rH 

to 

c*. 

to 

to 


Tf 

as 

OS 

LO 

LO 

LO 

LO 

aJ 

*h 

to 

to 

rH 

CM 

rH 

i-H 

rH 

o 

to 

CM 

CM 

CM 

CM 





rH 

1 

rH 

1 

rH 

1 

rH 

rH 

1 



1 

rH 

1 

rH 

1 

rH 

1 


fO 

SO 

LO 

i-H 

rH 

Os 

sO 

to 

to 

to 



00 

00 


so 

OS 

CM 

o 

o 

LO 

as 

LO 

sO 

LO 



CM 

LO 


rH 

l0 

00 

rH 

o 

CM 

00 

rH 

sO 




CM 

LO 

rH 

rH 

sO 

00 

rH 


to 

as 

to 

O 




rH 

sO 

CT3 

o 

O 

00 

LO 

o 

00 

c- 

00 

O 

CM 

L0 

OS 

rH 

O 


o 

o 

i-H 


00 

c- 

c- 


o 

o 

Ht 


00 

00 


00 

rH 

rH 

as 

as 

rH 

L0 

oo 

oo 

sO 

as 

f". 

to 

to 


to 

o 

LO 

LO 

00 

sO 

h~ 

rH 

t'' 

c- 

r-- 

to 

LO 

SO 


o 

rH 

as 

sO 

so 

to 

o 

to 

oo 

oo 

rH 

LO 


rH 

o 

o 

o 

CM 


sO 

CM 


LO 

o 

as 

LO 

sO 


© 


LO 

o 

o 

os 

rH 

CM 

CM 

CM 

o 

CM 

L0 

O 

rH 

rH 


o 

1 

rH 

1 

CM 

1 

CM 

1 

to 

1 

to 

t 

to 

to 

1 

o 

rH 

1 

CM 

1 

to 

1 

to 

1 

to 

1 










l''. 

i-H 

o 

LO 

rH 

LO 










00 

TT 

o 

US 

© 

as 










o 

OS 

o 

to 

to 

CM 

r-j 

CO 

LO 

o 

o 

o 

o 

o 

o 

o 

o 

CM 

L0 

LO 

to 

oo 

o 

f— 1 

CM 

to 

LO 

o 

LO 

o 

o 

rH 

CM 

to 









rH 

rH 

CM 








i 



40 


The chance of collision per unit time is equal to the probable number 
of molecules of the background gas whose centers lie within a cylinder of 
base area 4a and height V where a is the collision cross section of a particle, 

V the relative velocity. 

Let 6 c be the angle between c and c_' , and let <p be an azimuth angle for 

* 

The number of background molecules per unit volume for which c’, 9, 4> 
lie within small ranges dc' d9^ d<J> is 

. ~(c' 2 /vb . 

(n/Tr ' V )e c ' ^ sine d0 dcf>dc' 

m c c 

where v^ is the most probable thermal speed of the background gas molecules. 

Multiplying by 4aV and integrating over all values of <f> gives the number of 

background particles within the cylinder of volume 4aV such that c', 0 c lie 

within dc', d6 as 
c 

. -c c 2 /vb 

( 8no/v if )V e c' sine d6 d (2.11) 

m c c c 

when c, c' are given, V depends on 0 c as 

2 2 2 

V = c +c' -2cc' cosO 

c 

Differentiating with respect to 9 for constant c, c' 

VdV = cc' sin9 de 
c c 


Substituting in (2.11) gives 
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3 1/2 -C c,2 / V ^) 2 

(8no/vV /z )e m (c'/c)dc' V z 


dV . 


( 2 . 12 ) 


2 

Integrating V with respect to V, keeping c, c' constant 


c+C 

V 2 dV = V 3 /3|,^! 


c-c' 


2 2 2 
j c(c +3c' ) c' > c 

2 2 2 
J c ' (c' +3c ) c’ < c 


Thus integrating (2.12) with respect to V gives 
when 

c’>c ; (16 no 


3 1/2 - (C ' 2/V 5 2 2 

/3v 7T // )o m c'(c + 3c' ) dc' 


when 


c' <c 


_(c' 2 /v 2 ) 

(16 na/3vV /2 )e m (c ,2 /c)(c' 2 + 3c 2 ) dc’ 


The mean collision probability per unit time is then found by integrating 
from c' = 0 to ” using the appropriate expression as c' is greater or less 
than c 


16no 

_ 3 1/2 
3V TT 
m 


, 2 . 2 

2 2 _c /V 

c'(c +3c' )e m dc' 


c ,2, ,2 _ 2, 
c' (c 1 +3c ) 


-(c’W) 

m j i 

e dc f 


The first integral may be evaluated directly as 


(2cW m * 2/3) vV C /V ” 


The second integral cannot be evaluated in finite terms, but replacing 
c ,2 /v 2 by y 2 the integral becomes 
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v 5 f c / v m i 2 

_E y?(y 2 + 3c 2 /v )e y dy 
'o 


which after continued integrating by parts becomes 


v 5 f - (c 2 /V 2 ) 0 2 _ , , 2 

m I m' c ,2c A 3. A 3.2c , ... 

r r e V ( ~ 4-* 4 ( — + 1} 

v m v v 

m m 


c/v 2 
m e‘ y dy 


Hence the sum of the two integrals is 


3v 5 -c 2 /V 2 2 

m ,c ' m 2c .. 

5P- (- e * — + n 

m v 


c/v 

r m i. 


e ' y dy) 


This may be expressed as 


5 2 2 

3v° -c/v 1/2 0 2 

m ,c m tt ,2c , , .... 

4c~ v 6 + 2 ( — * 11 erf(c/v n,” 

m v 

m 


Where 


erf (C) = - 


172 


5 2 

-x , 
e dx 


the so called error function. 

Alternatively Jeans (1954) has tabulated values of the function 


-£ 2 2 -x 2 

<K€) = e ^ + (25 +1) e x dx 

} o 


(2.13) 


(2.14) 


In terms of which (2.13) becomes 
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3v 

-r-^- ^(c/v ) 
4c m 


The mean collision probability per unit time, v for a test particle of speed 
c is then given by 

v =* 4na (v^/ctt^^^)^(c/v ) . (2.15) 

The above derivation is for zero mass motion of the background gas. 

The Monte Carlo study requires a background gas of mean velocity v^, so that 
the above results may be used if c is replaced by an effective test velocity 
found by subtracting v^ from the true test velocity for the study, as 


c = 




‘ ^bi 


At this point, a quantity 

4v 2 

V rel ■ -nr (2 - 1 ' 

ir c 

is introduced, which may be regarded as the mean effective collision speed 
of the test and the background particles. 

It is assumed that the time fit representing any mo is so small that 
the probability of two collisions occurring therein is negligible, so that 
collisions with each background class are mutually exclusive. Hence total 
collision probab lity per unit time is found by addition and the total collision 
probability in time fit, P(bang) is 

P(bang) = (v 1 +v 2 +v 3 ) fit 
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f 


Subroutine VREL first determines the equivalent test particle speed 

I cl as above, and hence the current speed ratio S = |c|/v . If this is found 
— — m 

to be less than 0.1 the corresponding error function is taken equal to S, 
and if S is greater than 1.5 the error function is taken equal to 1. The 
expression (2.14) for ip (S) is then evaluated directly. Foi values of s between 
0.1 and 1.5, 15 values of the function ^(S) are stored for s from 0.1 to 1.5 
in steps of 0.1. The ip (S) corresponding to any S value is then used, and 
in all cases expression (2.16) is used to return the value of V re ^. 

2.4.5 Collision dynamics, subroutine BANG, PTNR3, VELS 

Having chosen the class of the colliding partner, it remains co determine 
the test molecule velocity after collision. 

One way to do this would be to select a collision partner from the 
partner class, according to some distribution. The. laws of dynamics could 
then be used to determine the test velocity after collision. 

However the distribution of partner molecules should be weighted by the 
relative velocity of the colliding molecules, so that the resulting expressions 
would be similar to those used in picking the initial velocity of a test 
molecule (Section 2.4.4). In this case the variable S z would be different 
for each collision. Since S z determines the coefficients used in the polynomial 
fit. a large table of possible coefficients would be needed. This approach 
would lead to very large penalities in time and storage. 

If the test molecule were completely accommodated into the partner class, 
its new velocity could be determined from the mean velocity of the partner 
class, plus a random velocity drawn from a Maxwellian distribution of thermal 
velocities corresponding to the partner class temperature. A relatively 
simple procedure. 
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However the test molecule cannot be regarded simply as a member of the 
partner class after collision. It is found that the original velocity tends 
to persist so that the expectation of the thermal velocity after collision 
is in the same direction as the original test velocity (Jeans, 1954). 

The test and partner molecules are assumed to be elastic spheres of 
cross section a. The directions of motion relative to the center of gravity 
of the two particles are AB, DE before impact, BC, EF after impact (Figure 2.9). 

In order for a collision to occur, AB produced must cut the plane perpen- 
dicular to AB within a circle of area 4o about E, say P. Also, all positions 
of P within the circle are equally probable so that the probability of EP 
being between r and r + dr is 


irrdr 

2a 


as 


1 

= sxn y 




1 . 

2 sm 


<J> dd> 
c c 


/■4ov 

r = ( r ) 


1/2 


sin j * c 


Thus all directions for EF are equally probable. Hence the expectation 
of the component of velocity of either molecule after impact in any direction 
is equal to the component of the mass center velocity in that direction. 

For molecules cf equal masses if OP and OQ represent the initial velocities 
of the test and colliding molecules, then OR represents the velocity of their 
mass center (Figure 2.10). 

To find the average velocity of the test molecule after impact OR is 
first averaged for all directions of the partner molecule, keeping its 


/ 



Figure 2.9 Collision model. 



48 


magnitude constant. The average component of OR perpendicular to OP is zero, 
leaving the average of ON to be determined. 

The probability of collision is proportional to the relative velocity 
of the test and partner molecules, thus the probability that POQ lies 
between 6 and 8+de is 


PQ sin 0 d6 


giving the average value of ON as 


fir 


ON = 


ON PQ sin6 d9 
fir 

PQ sin0 d0 
o 


let OP = c OQ = c' PQ = V 

so that 


V 2 = c 2 + c' 2 - 2cc' cos 0 


(2.17) 


Then 


ON = j (OP+OM) = j(c+c'cos0) 


1 2 ,2 ,, 2 . 

= 4 C (3c +c' -V ) 


differentiating (2.17) gives 


VdV = cc' sin0 d0 


OK - / . 2 -V 2 I V 2 dV 

4c / V dV 


giving 
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using the previous results (2.15) and integrating from V = |c'-c| to V * c*+c 


ON 


4 4 

15c + c» 

10c(3c 2 +c' 2 ) 


c > c' 


qm _ c(5c» 2 +3 c 2 ) 

2 2 
5 (3c* +c) 


c < c' 


The expectation of the velocity after collision is therefore in the 
same direction as the velocity before collision, for either molecule 


let c/c* = k and ON = a 


then a/c represents the persistence of original velocity in collision given 


by 


a/Cj^ 


4 

15k + 1 

10k 2 (3k 2 + 1) 


K > 1 


a/ C 


3< 2 + 5 
5(k 2 + 3) 


k < 1 


evaluation of these expressions for various values of k between 0 and °° 
show that a/c varies from .333 to .500. 

The distribution of k values is found to be 


5k (3tc +1) 

V + K J 


So that multiplying by the mean expression for a/c, and integrating for tc 
between 1 and « gives the mean persistence of all velocities after collision 
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°° 25k 4 + 6k 2 + 1 
1 4* / 2k(1 + < 2 ) Vl 


d< = I + i — lnCl+v^) = .406 
4 4/2 


The above derivation gives the mean persistence of the test molecule 
velocity after collision in the direction of the test molecule velocity before 
collision. The derivation was carried out assuming a background gas with zero 
mean mass motion. For the case of collisions with a mean motion of the back- 
ground gas, the system may be reduced to the above situation by adding to the 
test molecule velocity, a velocity equal in magnitude and opposite in direction 
to the mean velocity of the background class. 

The algorithm for generating velocities after collision must incorporate 
the concept of persistence of velocity, and must also be applicable in the 
following two extreme cases: 

(A) Test molecule at rest; warm partner gas with zero mean motion. 

(B) Test particle with finite velocity v , cold partner gas with zero 
mean motion. 

For case (A) the temperature of the test molecule population after 
collision may be found analytically as follows. 

Using the spherical coordinates defined in Section 2.4.4 the flux of 
molecules crossing the plane normal to the velocity of the collision partner 
is given by 


K 


■TT 

■2ir 

•o J 

0 


r - (v'Vv'b 2 

V cos 6 e m V 

o 


sin8 dV d<j> d9 


where K is a constant, V is the magnitude of a partner velocity, 0 is a polar 
angle with respect to the partner velocity vector and <f> is an azimuthal angle. 
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Integrating for 6 and 4 1 gives 


r 3 -cv 2 /v 2 ) 

K V e 2ir (1/2) dV 

o 


Hence the mean speed of partner molecules is given by 


r a 

1 V 4 e m dV 


f v 3 e "(V dV 


and the mean squared speed by 

„ 5 -(V 2 /v 2 ) 

-2 r 0 V 5 a m dV 

^ f v 3 e-^dV 
J 0 


2— j- = 2 v 2 

4 m 


Cl/2) v, 


m 


The mean squared speed of the gas at a point is given by 


/oo 


2 2 
4 -^2 


e dV , . 

o 3 

2 y 2q = 2 v m 


point f V 2 e ~ c yt/ 0 


m-'dV 


Thus 


V‘ = (4/3) V 2 . . 
flux v ' point 


The flux energy is divided equally, in the mean, at the first collision of 
type (A) . A test molecule initially at rest therefore acquires an average 
temperature 2/3 that of the background gas at a point in its first collision. 

In the case (B) the test molecule is scattered isotropically, and using 
the Jeans (1954) result for c/c ? = the persistence is 0.5. The velocity 
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of the test after collision is on a sphere of center 0.5 and radius 
0.5 iv I in velocity space, where is the original velocity of the test 
molecule. A warm partner gas in case (B) will lead to a thickening of the 
sphere to a spherical shell of thickness corresponding to a velocity drawn 
from a Maxwellian distribution of temperature equal to half that of the back- 
ground gas . 

The above considerations suggest the following method for generating 
the test molecule velocity after collision. The test velocity after collision 
v! is given by 


il - p(vV * 4 * (2 ' 18) 

where v^. , v^ are the test velocity and mean background class velocity 
before collision 

r is a random unit vector drawn from a spherically symmetric 
distribution 

[v^rT is a random thermal velocity drawn from a Maxwellian 

distribution with temperature rT, where T is the temperature 
of the background gas (see Section 2.4.6) 
p, q, r are constants whose values depend on the type of collision. 

The range of values of p, q, and r is as follows 
Case A Case B 

p 0.33 to 0.5 (the persistence) 

q ^ 0.5 for this problem 
r 0.67 to 0.5 


Energy conservation for the special case of collision between a test molecule 
of class 1 and a partner, also of class 1, gives an additional constraint; 
in this case 


4 = -1 


I t = Vj * [v r ] T 


where v^, T are the mean velocity and temperature of class 1 


Vj = P [v. r ] T + + + X 


•r J T' 


-r J rT 1 ^b 


= ptx r ] Ti - q[x r 3 Ti - [v r 3 rTi - ^ 

, 2 2 ’ 1/2 r , 

= (p + q + r) [v r ] T ^ + ^ 

But since the test molecules after collision must have the temperature of 


class 1 molecules, T^, 


(p 2 + q 2 + r) = 1 


For the purposes of the present study the values 


p = q = r = 0.5 

were chosen, as an adequate representation of all the above considerations. 

The subroutine PTNR3 uses subroutine VELS to give both [v ] T and r by 
independent random selections from a Maxwellian distribution. In the selection 
of r a software switch is set to 1, which causes VELS to jump round the 
selection of a velocity magnitude, so that only the direction unit vector 


is returned to PTNR3. 
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2.4.6 Thermal velocity with spherically symmetric distributions, 
subroutines VELS 

The fraction of particles moving outwards through unit area in unit 

3 

time (Figure 2.11) with velocities in d V about the velocity (V, 0, ip) is 


1 3 -( V /0 
- V e m simp dip de dV 

2ttv 

m 


(2.19) 


where v^ = most probable thermal speed of molecules. The marginal distribution 
functions f^, f ^ are found by integrating as follows (Fan, 1967) 


V* - 


f 2 » r 

o 


V 2 / 2 

’ 4 3 " V /v tn 

(l/2irv*) V e dv d9 sintp dip 


1 v m 


2ttv 


4 2 
m 


2tt simp dip 


= — simp dip 


f 0 de = 


,, 2 , 2 

r“ . - -V /v 

(1/2ttv ) V e m dV simp dip d0 
m 
o 


, v 4 

— 4 -f • 2 do 

2ttv 

m 


de 

2ir 


f v dv = 


2tt 


,.2 , 2 

fTT , , -V /V 

( 1 / 27 tv*) simp dip de V e dV 


-V 2 /v 2 

= (W2ttv 4 )V 3 e m dV 

,, 2 , 2 

/i , -V /v 

= (2/v 4 ) V 3 e m dV 
v nr 
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Again using equation 2.1, and setting R, uniformly distributed 
on the interval, 0 to 1, equal to the cumulative distribution function 
(see Section 2.1) 


rx 


R. 

1 


fCO 


where the random variable x has the distribution f(x) .which is then solved 
for x given any value of R^ . 

The above distributions give, for 6, V 


R e -17 

R t * 7 (1 - cosy) 

*V ■ <1 * e 


where R Q , R , Ry are independent random numbers uniformly distributed on 0 to 
1. The latter expression cannot be rearranged to yield V as an explicit 
function of Ry, however if the distribution function is re-stated in cylindrical 
coordinates, an expression may be obtained for V explicitly in terms of two 
independent random numbers Ry., , Ry r uniformly distributed on the interval 
0 to 1 (Pertmutter, 1966) 
whence 

V - v m nn < 1 ' ,R v ! R v r ) ‘ 2 - 20 > 

■ cos 2 (1 - 2 R ) 

6 = 2ir R 

V 

as used by subroutine VELS. 
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2.4.7 Surface interaction subroutine HITSPH 

As remarked above (section 2.3.5) a specular reflection model is assumed 
at the surface of the probe. Other workers (Fan 1S67) have used a more real- 
istic diffuse reflection model incorporating Lambert’s cosine law of scattering. 
Here the velocity and temperature of the re-emitted molecule is determined by 
its incident velocity, probe geometry and probe temperature. It is felt that 
this additional refinement can be traded against the extra storage required 
for variable class 2 temperatures, without severely compromising the realism of 
the simulation. 

For specular reflection the molecule temperature is assumed to be unchanged 
by reflection. The test molecule velocity after reflection is found by re- 
versing the velocity component radial to the sphere, that is the velocity after 
reflection is given by 


V 2 ( Vi)i (2 - 21) 

where is the test velocity before reflection 

r is a unit vector radial to the sphere at the point of impact. 

2.5 Time-Saving Techniques 

In any program involving multiple iterations, the problem of time economy 
becomes very important. The Monte Carlo program obviously falls into this cate- 
gory, since in tracing one particle, parts of the program are executed several 
hundred times, and a typical iteration consists of 10,000 particles. 

An examination of execution times for various FORTRAN operations on the 
PDP-15 points to some obvious economics. 

Real exponentiation on the PDP-15 takes 30 times as long as real multi- 
plication, so that raising to any integer power less than 30 should be accomplished 
bv reneated multinlication . rather than exoonentiation. 
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The function SQRT is four times faster than real exponentiation, so that 
it should be used in preference to exponentiation to power U.5. Logical IF 
statements are faster in execution than arithmetic I! statements, so that where 
only two branches are desired a logical IF should always be used. Complex logical 
IF statements such as 


IF(A.GT.B.OR.A.GT.C) GO TO 1 

should be replaced by sequences of simple IFs as 

IF(A.GT.B) GO TO I 

IF(A.GT.C) GO TO 1 

as the complex expression will always be tested in its entirety whereas the 
first simple expression may fulfill the condition without the second being 
executed, especially if the condition most likely to be fulfilled is placed first. 
An expression such a. 


2 3 4 

y = a,x+a„x +a,x +a.x 
12 3 4 

should be arranged 

y = (a 1 +[a 2 +(a 3 +a 4 x)x]x)x 

The former requires ten multiplications and three additions, while the latter 
requires only four multiplications and three additions. 

If the same array clement is to be used more than orce, it is quicker to 
set a simple (unsubscripted) variable equal to the subscripted array element, 
and use the simple variable in subsequent expressions, so that the array must 
be accessed only once. 
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Similarly if the same combination of constants is to be used more than 
once, a new constant should be introduced, equal to the combination. The above 
are good practice in any FORTRAN program, although the time saving will be 
greatest where much of the time is spent in execution, rather than in input/ 
output . 

2.6 Operating the Program 

The following handler assignment must be made before calling the loader; 

"A DTC2 -4,-5/DTCO -1/DTF3 3" 

The C handler is limited to read only (for program loading) and is much smaller 
than the usual A handler (680^ instead of 2290^). The F handler is also small 
and will read or write in the non-file oriented mode (for data). The loader is 
now called by typing "GLOAD" and when the loader responds, the main programs are 
loaded by typing "+MCSPH7". All subroutines are contained in a .LIBR5 BIN user 
library file, and are loaded automatically. 

The main programs request the input data shown in Table 2.2 which must be 
determined by the user, and typed in on the teletype. 

The number density, temperature and velocity for free stream are determined 
by the problem to be simulated. The number of cells is determined by core avail- 
ability, and IZMAX and IRMAX have been set for the present program to maximum 
figures of 30 and ten respectively. 

The cell sizes must be determined by short trial runs. A good criterion is 
that the cell size should be of the order of half a mean free path, and the time 
for one mo such that the distance travelled in a mo is half an axial cell, based 
or the mean effective speed in the free stream. The position of the probe 
(ZPROBE) is constrained by the fact that class 2 velocities are stored only for 
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Symbol 

KM 

IZMAX 

IRMAX 

CREF 

DU 

DZ 

DR 

DT 

ZPROBE 

N 

T1 

VI 


FMFP 


T1S 


FMFP 

T1S 

LI 

IRT 

NI 


TABLE 2.2 Input data for programs 
Main Programs 

Meaning 

Number of particles per iteration 
Number of axial cells (30) 

Number of radial cells (10) 

Reflection coefficient (normally = 1) 

Seed for RANDU, see section 2.4.1 
Axial cell size 
Radial cell size 
Time increment per mo 
z coordinate of probe front 
free stream number density 
free stream temperature 
free stream effective velocity 

Data reduction programs, NEWBG7 

free stream mean free path 

free stream stagnation temperature 

Standard deviation program STDV7 

free stream mean free path 

free stream stagnation temperature 

switch; 0, 1, 2 or 3 for density, temperature, 
axial or radial velocities respectively 

radial cell number 

iteration number 
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the front ten axial cells from the exit plane. The front tip of the probe should 
not extend more than seven axial cells from the exit plane. 

The reflection coefficient CREF is always set 1 in this form of the program. 

The number of particles per iteration is set by trial and error after per- 
forming standard deviation checks (see section 2.3.10). A sample size of 10,000 
is adequate in most cases. 

After typing in the data, the program will repeat all parameters as a check. 

If the data is found to be incorrect, striking CTRL/S will restart the program and 
the whole must be retyped. Otherwise execution will begin. The program will con- 
tinue until it is dumped or control is returned to the monitor by striking CTRL/C 
on the teletype. 

After an adequate number of iterations have been completed and stored on 
tape, the data may be retrieved and printed by loading the off-line data re- 
duction program NEWBG7, and either RESNM7 or RESRW7. In either case, the F 
handler must first be assigned by typing "A DTF3 3". The loader is then called 
by typing "GLOAD" and the programs loaded by typing "-‘-NEWBG?, RESNM7" or 
'W-NEWBG7, RESRW7 . 

Alternatively, the programs STDV7 and RSTDV7 may be used to calculate stand- 
ard deviations for any parameters by typing: "A DTF3 3 M , "GLOAD", "«-STDV7, RSTDV7". 

The off-line data reduction programs also require input parameters from the 
teletype (Table 2.2). FMFP and T1S are the free stream mean free path and stag- 
nation temperature used to normalize output parameters. LI, IRT, NI determine 
which standard deviations will be computed. 
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3. RESULTS 


3. 1 Convergence and Accuracy 

For any iterative scheme, a criterion must be found to evaluate the conver- 
gence rate and determine a satisfactory end point of the procedure. In the present 
study a continuous plot was made for the density ratio in a cell in the flow field 
after each iteration. Figures 3.1 and 3.2 show that the density ratios for one 
particular cell (iz =23, ir = 5) at each altitude, converge satisfactorily to 
stable values. After the flow field has become stable there remain statistical 
fluctuations, so that error limits must be established such that the iterative 
procedure may be stopped when further changes between successive iterations are 
lost in statistical fluctuations. The Monte Carlo scheme permits an easy solution 
of this problem, as the standard deviation of any flow-field parameter may be 
estimated in the following way. An iteration of, say 10,000 molecules is com- 
pleted in the usual way and the results are stored on tape. 

The same iteration is then repeated, say ten times, using short runs of 
a lesser number, n of molecules, say 1000, the resulting parameters from each 
short run are also stored and from these, separate programs (STDV7, RESNM7) compute 
a best estimate of the standard deviation of any required flow field parameter, 
X, according to the usual formula 


= { EX 2 /n- (IX/n) 2 }n/n- 1 


(3.1) 


It may be shown that the standard deviation of a Monte Carlo calculation 


decreases as 1/t/jT" where N q is the sample size, so that a best estimate of the 


lOn 


for the original run of lOn molecules will be 


a = o //To" 

lOn n 


standard deviation a 
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The best estimate for the mean, )T, of the lOn run will be the same as 
that of the n run. 



Wherever possible the plots that follow are given error bars corresponding to 
plus and minus one standard deviation computed as above. 

Since for a normal distribution, 64% of values fall within plus and minus 
one standard deviation of the mean, changes of the order one to two standard 
deviations between successive iterations are considered statistical fluctuations. 
Figure 3.1 shows the mean and standard deviation computed for the seventh iter- 
ation for a sample size 10,000 and indicates that for these conditions corres- 
ponding to 70 km, five iterations would be adequate. Figure 3.2 shows a similar 
plot for 75 km where 12 iterations are required. As ambient density increases, 
the effect of class 3/class 3 collisions becomes more important, so that more 
iterations are required to establish the final flow field. For 75 km, the first 
13 iterations had a sample size of 2000. since statistical fluctuations would at 
first be lost in flow field convergence. Iteration 14 was of 10,000 molecules, 
to give low statistical fluctuation in the final results. 

3.2 Flow Field 

Figures 3.3 and 3.4 show density contours for the two altitudes treated, 70 
and 75 km. Points on the contours are interpolated from plots of density ratio, 

P/ p co a g a i nst radial distance tor each axial cell. Near the axis, accuracy of the 
data is poor because of the small cell area in the radial plane, and correspondingly 
small sample size. A quadratic curve fit is used for density ratio as a function 
of radius since its Taylor expansion reduces by symmetry arguments, to that of an 
even function, having only even power terms. Density ratio in the first three 
radial cells was plotted against r and a best straight line drawn through the 
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resulting points. Figure 3.5 shows density contours for a sphere at K n = 0.26, 

S = 5 from Bird (1968), a Monte Carlo simulation. The form of the contours is 
similar to those of the present study, Figure 3.4, but with the shock front 
spreading out further from the probe axis and swept back more steeply because of 
the higher speed rates, S. 

A quadratic fit was used also to give points on the front stagnation line 
density profiles shown in Figures 3.6 and 3.7. Error bars are not shown on these 
two plots. Figure 3.6 also shows the results of electron beam density measure- 
ments for a sphere taken under conditions similar, though not identical, to those 
of the Monte Carlo run. The measurements were made by Russell of JPL and are 
taken from a paper by Vogenitz ejt al . (1968). The experimental results show a 
thicker shock layer and higher density rise which would be expected for the lower 
Knudsen number and higher speed ratio used. Exactly corresponding experimental 
results are not available at the present time. 

Figures 3.8 and 3.9 show velocity vector plots in the plane of the sphere 
axis. Stagnation near the head-on point of the sphere, and deflection of the 
flow off axis may be clearly seen. Figures 3.10 and 3.11 show stream-lines 
sketched from these velocity vectors which show the expected form for specular 
ref lection, approaching a tangent to the sphere surface. 

Figures 3.12 and 3.13 show plots of the normalized mean temperature ratio, 
(T-TJ/CT -T w ), along the stagnation line at the two altitudes. Where T = mean 
temperature computed ‘rom the average energy over all classes relative to an ob- 
server moving with the mean velocity of the flow in the cell 

T^ = static temperature in free stream 

2 

T goo = stagnation temperature in free stream given by T Soo = T^ (l+l/2(y-l)M tE ] 
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Figure 3.fc 70 km stagnation line density profile, 
comparison with experiment. 
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Figure 3.12 70 km stagnation line temperature profile 



Figure 3.13 




where y = 1.66 the ratio of specific heats for the monatomic molecules of the 
simulation 

M = the free stream Mach number 

00 

Comparison of Figures 3.6 and 3.12, 3.7 and 3.13 show that a sharp tempera- 
ture rise coincides with the beginning of the shock, represented by the density 
rise . 

Figures 3.14 and 3.15 show density ratio profiles off axis, plotted directly 
from the data for the fifth radial cell. Error bars are shown, representing plus 
and minus one standard deviation calculated as in section 3.1. The fifth radial 
cell represents the limit of the sphere surface at both heights, so that the 
fall of density behind the sphere may be seen. Large standard deviations near 
the surface correspond to cells cut by the surface, and thus having less samples. 
Densities have been corrected to allow for the reduced cell volume corresponding 
to the cut cells . 

Figures 3.16 and 3.17 show off axis plots for 70 km of temperature ratios 
as previously defined, and axial and radial velocity ratio, again for the fifth 
radial cell. The temperature in Figure 3.16 rises at the shock, then falls back 
towards free-stream conditions behind the snhere, an effect somewhat blurred by 
high sampling errors in the cells cut by the sphere. Figure 3.17 shows the fall 
of axial velocity and corresponding rise in radial velocity as the particles are 
deflected past the sphere, again returning towards free-stream conditions behind 
the sphere. 

3 . 3 Drag Coefficient 

The sphere drag coefficient was calculated for the two altitudes 70 and 
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where AP = change of momentum flux ac oss the sphere 
o - free stream dens i tv 

oo 

V = free scream velocity 
A = projected frontal area of sphere 
AP may be calculated from the density ratio and velocity 
Monte Carlo program: 


10 2 ■> 

AP = E in { (N.V. A.). -(N.V. "A.) 1 

j =1 J J J tn j j j out 


ratio given by the 


where m is the mass of a molecule, A is the cell area in the radial plane, N 
and V are the number density and mean velocity and subscript j refers to the 
jth radial cell 


(3.2) 


where 


C n = (2 dr/r , )[E. -2 J 

D sphere L m out J 


10 • v- 2 

1 = 1 £1 (11) (2 j - 1) 

. , P 00 Voo j 
1=1 

r sphere = s P“ ers radius 

dr = radial cell width 

The resulting values of C D are shown in Table 3.1, along with values from 
Aroesty (1963) . 

The Reynolds number Re. the usual parameter for experimental work is defined 


as 

pa>Voo2r , 

R P = *P he - re 

Uoo 

where Uoo = free stream coefficient of vicosity. 
number, M, and Knudsen number, Kn by 


Re = 1.37 


Kn 


Re is related to the Mach 


Drag coefficients would be expected to be higher than experimental values; 
because the calculation makes no allowance for the axial momentum flux lost through 
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TABLE 3.1 Comparison of experimental and calculated 
drag coefficients. 

Parameter Height (km) 


70 73 

Kn 0.1 0.21 

M 2.74 2 . 74 

Re 37.6 17.9 

C D (Monte Carlo) 4.8 4.2 

C D (Aroesty 1963) 1.82 2.23 
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the cylindrical walls of the system with class 
pectation is confirmed by the results in Table 


2 and 3 particles . 
3.1. 


This ex- 
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4. CONCLUSIONS 

4 . 1 Limitations and Suggested Refinements 

The major simplifications involved in the present study relate to the 
mechanism of intermolecular collision and surface interaction. Hard sphere 
molecules have been used in the present program, for the sake of simplicity, 
but the program could easily be modified to model other types of intermclecular 
force law with different exponents v. Hard sphere molecules correspond to 
v = °°, while for example Maxwellian molecules correspond to v = 5 and another pro- 
posed model has v = 9 (Bird 1970). The principal difference between these models 
relates to the shock thickness for very strong shock waves, which is predicted 
to be independent of shock strength for hard spheres and directly proportional 
to shock Mach number for Maxwellian molecules. 

The present study treats only monatomic molecules, so that an obvious refine- 
ment would be the introduction of more realistic diatomic molecules, with their 
associated rotational and vibrational energy states. Again this would not change 
the program in overall concept, although execution and storage would be increased 
somewhat . 

Specular reflection at the probe surface, while simple to program, is not 
very close to known surface interactions in real studies, A more realistic model 
involves diffuse re-emission of the incident molecules, with velocities dependent 
on surface temperature as well as incident velocity. 

Further refinement of the program using core storage at present available on 
the PDP-15 will probably involve splitting the cell system so that only part of 
the cells are resident in core at any time, the remainder being stored on DECTAPE. 
The system might be divided into two cylindrical shells, each containing five 
radial cells. Background parameters for the inner block would be read into core 
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first, and molecules would be introduced at the entrance plane. They would 
be traced in the usual way, except that on crossing the radial boundary into 
the outer block, their positions and velocities would be stored for use in 
the second part of the iteration. After a suitable number of tests had been 
run, the accumulated data for the inner cells would be read onto tape, and 
the background parameters for the outer cells read into core; "swapping". 

Molecules would now be introduced, from the input plane, and a suitable number 
recorded, including those crossing back to the inner block. The molecules 
which had previously crossed to the outer block would now be re-started at their 
stored positions and velocities. Most of them should leave the end wall of the 
system, while a lesser number would be stored as returning to the inner block. 

The blocks would now be swapped again and the new molecules, stored at the boun- 
dary, released. After two or three swaps the number of molecules to be intro- 
duced at the dividing boundary should have fallen to an insignificant level, 
when the iteration would be completed in the usual way. 

Because of the small area of the innermost cells, few molecules are introduced 
there, causing higher statistical fluctuations. A system of weightings by radius 
could be introduced to improve accuracy of the near stagnation line parameters, 
while maintaining the effectively uniform distribution with area. 

4 . 2 Summary 

The Monte Carlo direct simulation technique has been used to model supersonic 
(M = 2.7) gas flow round a spherical probe in the transition regime. The two 
cases completed correspond to a probe of diameter 1 cm at heights of 70 and 75 km, 
the D region of the ionosphere. Monatomic hard sphere molecules are used and 
specular reflection is assumed at the probe surface. The technique can be modi- 
fied to treat other types of molecule and surface interaction without undue dif- 


fi culty . 
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Density contours are shown for the two altitudes, and compared with those 
of an independent Monte Carlo simulation. Velocity vector fields and stagna- 
tion line profiles of density, temperature and velocity are compared with avail- 
able experimental data. Drag coefficients are computed and compared with those 
derived from wind tunnel testing under similar conditions. Standard deviation 
estimates are included for data plotted directly from the program output. 

The work was carried out on the Aeronomy Laboratory PDP-15 computer, a 
small (16K) machine used principally for on-line data processing of a partial- 
reflection experiment. The study shows that Monte Carlo simulations are feasible 
using the limited core memory, often available in on-line machines, having a 
low duty cycle. 

Because of the nature of partial-reflection measurements, the PDP-15 on- 
line processing takes place during daylight hours only, so that the machine is 
available at night for Monte Carlo runs. A typical iteration of 10,000 mole- 
cules takes five hours for an altitude of 75 km. The results of each iteration 
are stored on DECTAPE. The present Monte Carlo results will be incorporated 
in future work to simulate ion-collection by the rocket-borne probes launched 
at Wallops Island by the Aeronomy Laboratory. These simulations will relate 
probe currents to ambient conditions in the D region, a problem intractable using 
present analytical methods. 
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APPENDIX 


MCSPH7 

A DTC2 -4,- 5/DTC0 -1/DTF3 3 

MONTE CARLO SIMULATION OF GAS FLOW PAST A SPHERE IN TRANSITION 
REGIME. 


(M.X.S. 

UNITS) 

KM: 

# PARTICLES /I TER ATI ON 

IZMAX: 

#AXIAL CELLS 

IRMAX: 

FRADIAL CELLS 

CREF: 

REFLECTION COEFFICIENT (1.) 

DU: 

SEED FOR RANDU (22g7.) 

DZ: 

AXIAL CELL SIZE 

DR: 

RADIAL CELL SIZE 

DT: 

TIME INCREME NT /MO 

ZPROBE: 

AXIAL DISTANCE FROM ENTRANCE PLANE TO FRONT OF PROBE 
(PR OBE RADIUS:. 005) 

N: 

FREE STREAM NUMBER DENSITY 

Tl: 

" " temperature 

VI: 

• " VELOCITY 


DATA SWITCH SET EFFECT 


00 

WRITE TEST CLASS & CELL 

P EACH MO 

01 

STANDARD DEVIATION MODE 

AFTER THIS ITERATION 

05 

PREPARE TO DUMP 


06 

RESTART WITH BACKGROUND 

FROM DECTAPE 

14 

REPEAT ITERATION 1 


17 

WRITE DATA ON DT3 


I NTEGER 

CT,CP,VZ3,VR3,VZ2,VR2 


REAL N 




LOGICAL LSSW 

DIMENSION N 1 (30, 10), N2 (30, 13), N3 (30, 10) ,M 1(30, 10) ,M2(30,!0) 
2,M3(30,!0),SVZ3(20,10) t VZ3(23,10),SVZ3S(20,10>, 
3SVXY3S(20,10),IT3(20,13),YI (3), 

4SVR3 (20, 13). VR3(20, 13), M2 0(10), M3 0(10) 

COMMON/II/RMAX,PI,IRMAX,VM/JJ/BT/MM/IZMAX,ZMAX/HN/M3,SVXY3S, 
2SVZ3S,SVZ3,SVR3,IZ,IR,N1 ,N2 ,N3 ,M1 ,M2 ,VZ3 ,VR3,IT3 , 

3 M10(10),FINT, FPEAL/KK/ZCTR ,R PR03E/LL/SVR2 (10,13), 

4 SVZ2(13,IO),VZ2(10,13),VR2(10,10)/III/RK2,VMT1 ,UZ2,UR2 

5 ,UZ,UR,V1 ,C0STH,T3C,SI NTH/JJJ/DR ,DZ,T 1 ,1 ,K,ZPROBE 

6 /KXK/IZC,CT,CP,S, N/LLL/L1 /KKM/KM/V/VZ,VX,VY 
CNU(ND,VA!,VA2,VA3,VBI ,VB2,C0STH,SINTH)=S*FREAL*FL0AT(ND)*VREL( 
2 VA l ,VA2 ,VA3 ,VB 1 ,V32,C0STH,SI NTH) 

WRITE (6,2) 

FORMAT (2X, 64H TYPE KM,IZMAX,IRMAX,CREF,DU,DZ,DR,DT,ZPPOBE 
2,N,Tt,Vl 1 PER LINE ) 

READ (4,4)KM,IZMAX,IRMAX,CFEF,DU,DZ,DR,DELT,ZPR0BE,N,TI,V! 
FORMAT(I6/I3/I3/F6.2/F6.3/(EI0.2)) 

WRITE (6,4) KM,IZMAX,IRMAX,CREF,DU,DZ,DR,DELT,ZPR0BE,N,T1 ,V1 
RsPAND'J(DU) 

I =0 

XMT=KM/10 
L 1:0 

REWIND 3 
CALL SETUP(N) 

IF(X.GE.KM) GO TO 416 
IF(IT.GE.10) L1=0 


NEW PARTICLE 
10 CALL NEWPOS'Z.R, THETA) 

COSTH=COS(THETA> 
SINTH=SIN(THETA) 
X=P*COSTH 
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Y:R*SlNTH 

IZ=Z/DZ+1. 

IZC=IZMAX-IZ+1 

IR=R/DR+I. 

MI0(IR):MI0(IR)+1 

CT:l 

VMrVMTl 

PALL «El'V(V| ,0.,VX,VY,VZ> 

DT=DELT*RA NDU (DU ) 

C 

C I NCRFKENT TIME! 

?l CALL INCPOS(Z,X,Y,R,COSTH,SINTH,VZ,VX,VY> 

DTsDELT 

C tlEXT CELL? 

IZNEW=Z/DZ+I . 

IR NEW~R /DR+I , 

IF(LSSW(0))WR*ITE(6,50)CT,IZNEW,IRNEV 
50 FORMAT (313) 

IFdZNEW, NE.IZ) GO TO 100 
IF(lRVEW.fJE.lR) GO TO 100 
IF(IZ.EQ.l) GO TO ! 00 
GO TO UR 
C MEW CELL 
100 IZiIZUEW 

IR=IRNEW 

CALL NEWCEL(Z,X,Y,R ,VX,VY,VZ,CREF,J) 

GO TO (300,110,120), J 
110 VM = VMT1 

C 

C COLLISION PROBABILITIES 

CNU1=CNU(NI(IZ,IR> ,VZ,VX,VY,V1,0.,1.,0.) 
IF(N2(IZ,IR).GE,1) GO TO 113 
IF(N3(IZ,IR).GE.l) GO TO 113 
PBA NG=C NU 1 *DT 
CNU2 = 0. 

CNU3 = 0. 

00 TO 118 

113 CNU2=CNU(N2(IZ,IR),VZ,VX,VY,UZ2,UR2,C0STH,SINTH> 

VM:S0RT(RK2*T3C) 

CNU3=CNU(N3(IZ,IR),VZ,VX,VY,UZ,UR ,C0STH,SI NTH) 
PBA NG= (CNU 1 +CNU2+CNU3 )*DT 
118 IF(RANDU (DU) .LT.PBA MG )G0 TO 140 

120 IF(MI (IZ,IR).GT. 130000) GO TO 312 

IF(M2(IZ,IR>.GT. 130000) GO TO 312 

IF(M3(IZ,IR).GT. 130300) GO TO 312 

130 CALL 3KXEEP(VX,VY, VZ) 

60 TO 21 
C 

C COLLISION: WHAT CLASS? 

140 IF (CNU2 .GT.0.) GO TO 150 

1 F(C NU3 «GT, 0, ) GO TO 150 
CP = 1 

GO TO 152 

150 YI(I>=CNUt*DT/P3AN0 

YI(2)=CNU2*DT/P9ANG 
YI(3)=CNU3*DT/PBANG 
CP = MONTE(YI) 

IF(CT.EQ.3) GO TO 160 
152 IF(CP.EQ.CT) GO TO 120 

160 CALL BANG(CT.CP) 

GO TO 130 

300 I F (LSSW (5 ) ) CALL DUMP 

310 X sK+1 

C 

C ENOUGH PARTICLES? 

IF (K.GE.KM) GO TO 312 
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IF(L I ,E0 • l ) GO TO 500 
GO TO 10 
C 

C NEW BACKGROUND; OUTPUT 
312 1=1+1 

VRITE<6,415> I ,IT ,K 

413 FORMAT (3H I = , 13 , 1 H . ,I2/3H Ks.IS) 

CALL RESULT <IZMAX,I RMAX) 

C LI SETS TO STANDARD DEVIATION MODE 
L! = 0 

IF <LSSW < 1 )) L I = 1 

416 CALL NEWBG(IZMAX,IRMAX> 

IT=0 

IF (LSSU( 14)) GO TO 9 

417 K =0 

GO TO 10 

500 IF(K.LT.KMT) GO TO 10 

CALL RESULT<IZMAX ? IRMAX> 

IT=IT+1 

WRITE CC|A15) I » I T ,K 
K = 0 

GO TO 9 
END 


C BANG7 

C SUPPLIES CORRECT COMBINATION OF PARAMETERS TO PTNR3 
C 

SUBROUTINE BANG (CT , CP) 

INTEGER CT,CP 

C OMMON/II I/RK2 , VMT 1,UZ2,UR2,UZ,UR,V1 f COSTH ,T3C,SI NTH/II /RMAX.PI 
2 IRMAX, VM/V/VZ , VX,VY 

C SET VEL FOR PTNR3;CP=CT DEALT WITH ABOVE UNLESS CT=3 
IF(CT-2)155 f 175,195 
J55 IF (CP ,GT,2) GO TO 170 

VM=.707*VMT1 

CALL PTNR3 (UZ2 ,UR2,C0STH,SI NTH,V1,0 # ,1,,0.) 

GO TO 210 

170 VM=.707*SQRT(RK2*T3C> 

CALL PTNR3 (UZ ,UR ,COSTH,$I NTH,V1 ,0., I , , 0, > 

GO TO 210 

175 IF(CP.GT.l) GO TO 1*5 

VM=.707*VMT1 

CALL PTMR3(V1,0.,1.,0.,UZ2,UR2,COS7H,SINTH> 

GO TO 210 

185 VMr.707*SQRT(RK2*T3C> 

CALL PTNR3 (UZ ,UR »COSTH,SI NTH,UZ2,UR2,C0STH,SI NTH) 

GO TO 210 

195 I F (CP-2) 197,199,205 

197 VMs »707*VMT I 

CALL PTNR3CV1,0.,1.,0.,UZ,UR,COSTH,SIHTH) 

GO TO 210 

199 VM=.707*WTI 

CALL PTNR3 (UZ2 ,UR2 ,C OSTH ,SI NTH ,UZ ,UR ,COSTH ,SI NTH ) 

GO TO 210 

205 VMr.707*SQRT(RK2*T3C> 

CALL PTNR3 (UZ ,UR ,COSTH,SI NTH,UZ , UR, C OSTH, SI NTH) 

210 CT=3 

RETURN 

END 



C BKKP7 

C ACCUMULATES PARAMETERS 
C 

SUBROUTINE BKKEEP (VX,VY,VZ) 

INTEGER CT ,CP, VZ3,VR3,VZ2,VR2 
REAL N 
LOGICAL LSSW 

DIMENSION Nl(30,l0),N2(30,l0),N3(30,10>,Ml (30,!0>,M2(30,!0) 
2, M3 (30,10) ,SVZ3 (20, !0) ,VZ3 (20,1 0) ,SVZ3S(20, 1 0) , 

3SVXY3S(20, 1 0) ,IT3(20,10),YI(3), 

4SVR3 (20,1 0) , VR3 (20, 1 0) ,M20( 1 0) ,M3 0( 1 0) 

C0MM0N/II/RMAX,PI,IRMAX,VM/JJ/DT/MM/IZMAX,ZMAX/NN/M3,SVXY3S, 
2SVZ3S,SVZ3,SVR3,IZ,IR,N!,N2,N3 ,MI ,M2,VZ3,VR3,IT3, 

3 M10(10),FI NT.FREAL/KK/ZCTR ,RPR0BE/LL/SVR2 (10,10), 

4 SVZ2(10,10>,VZ2(I0,!0),VR2(t0,10)/III/RX2,VMTl ,UZ2,UR2 

5 ,UZ,UR,V1,C0STH,T3C,SINTH/JJJ/DR,DZ,T!,! -K.ZPROBE 

6 /XKX/IZC,CT,CP 
IF(2-CT) 133,125,122 

122 MI(IZ,IR)=Ml(IZ,IR)+l 

RETUR N 

125 M2(IZ,IR)=M2(IZ,IR)+1 

IF (IZC*GT. 1 0) RETURN 

SVR2(IZC,IR)=SVR2(IZC,IR)+VX*C0STH+VY*S1NTH 

SVZ2(IZC,IR)=SVZ2(IZC,IR)+VZ 

RETUR N 

130 M3(IZ,IR)=M3(IZ,IR)+I 

IF (IZC.GT.20) RETURN 
VR = VX*C 0STH+VY*S1 NTH 
SVXY3S(IZC,IR)=SVXY3S(IZC,IR)+VR*VR 
SVZ3S(IZC,IR)=SVZ3S(IZC,IR)+VZ*VZ 
SVZ3(IZC,IR)=SVZ3(IZC,IR)+VZ 
SVR3(IZC,IR)=SVR3 (IZC,IR)+VR 
RETURN 
END 


C DMPTK5 
C DUMP S, RECOVER 
C 

SUBROUTINE DUMP 

COMMON/.' 'J/DR,DZ,T1,I,K,-PR0BE/KKK/IZC,CT,CP,S,N 
VR I TE > o , 5 ) 

5 FORMA'. ,£'A mount TAPE & TYPE* TQ2) 

CALL PAWSE 

C REPOSITION DATA TAPE 
REWIND 3 
DO 8 IT=!,I 
DO g IREAD=I,I6 
8 READ (3) A 

VRITE(G,10) 

10 FORMAT (22H PROGRAM HAS RESTARTED) 

RETURN 

END 


r HTTCDlC 

C SURFACE REFLECTION 
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SUBROUTINE HITSPH(VX,VY,VZ,Z ,R ,X,Y> 
COMMON/KX/ZCTR ,RPROBE/MM/IZMAX,ZMAX 
C ASSUME HIT AT PRESENT R 

ZSP=SQRT(RPROBE*RPROBE-R*R) 

ZSP=SIGN(ZSP, (Z-ZCTR) ) 

I F( (ZCTF+ZSP) .GE.ZMAX) ZSP=-ZSP 
Z:ZCTR+ZSP 

RM=SORT (X*X+Y*Y+ZSP*ZSP> 

C UNIT VECTOR FROM PR(*E CENTER TO POINT OF IMPACT 
RXU=X/RM 
RYU=Y/RM 
RZU=ZSP/RM 

VD OTR U = VX*R XU+ VY*R YU+ VZ*RZ U 
VX = VX-2 .* VD OTR U* R XU 
, VY=VY-2.*VD0TRU*RYU 
VZ=VZ-2.*VD0TRU*RZl' 

RETURN 

END 


| 

I 


C GIVF5 

C ACCUMULATED PARAMETERS TO DZCTAPE 
C 

SUBROUTINE GI VE (DR ,DZ , N) 

INTEGER VR2,VZ2,VR3,VZ3 
REAL N 
LOGICAL LSSW 

COMMON/I I /RMAX,P I »IRMA X,VM /MM/IZMA X,ZMA X/NN/M3 (30,10), 

2 SVXY3S(20,10),SVZ3S(23,13),SVZ3(20,10),SVR3(20, 10), 

3 IZ , IR , N1 (33,10),N2(30,10),N3(30,13),M1 (30,13), 

4 M2(30,10),VZ3(20,13),VR3(23,10),IT3(20,I0), 

5 Ml 0( 13) ,FI NT ,FREAL/LL /SVR2 (13,10), 

6 SVZ2 (13,10) ,VZ2 (13,13) ,VR2 (I0,13)/KK/ZCTR ,R PROBE /I II /RK2, 

7 VMT! ,UZ2 ,UR2 ,UZ ,UR ,V1 , THETA ,T3C/JJJ/DRD ,DZD ,T 1 
WRITE(6, 1 0) 

10 FORMAT (2BH PREPARE DT3 , THEN SWITCH 17) 

11 CONTINUE 

I F ( , NOT .LSSW ( 1 7) ) GO TO 11 

J = 1 

K = I5 

K2| = l 

K22:5 

K3 1 = 1 

K32=10 

C FOR EFFICIENT USE OF DECTAPE 
DO 20 L=I,2 

WRITE'3)((MI (IZ,IR),IZ:J,K),IR=1,10),IZMAX 
WRITE «.3) ( (M2 (IZ ,1 R ) ,1 ZrJ ,K) ,IR=1,I0> 

WRITE (3) ( (M3 (IZ,IR),IZ=J,X> ,IRr 1 ,10),IRMAX 
WRITE(3)((SVR2(IZ,IR),SVZ2(IZ,IR),IZ:K21,K22),IR:1,10),PI 
WRITE (3) ( (SVR3 (IZ, IR ) ,IZ:K31,K32),1R= 1,10), R MAX 
WRITE(3)((SVZ3(IZ,IR),IZ=K3I,K32),IR=1,10) 
'•RITE(3)((SVZ3S(IZ,IR),IZ = K31 ,K32),IR=1,1 0),ZMAX 
WRITE (3) ( (SVXY3S(IZ,IR) ,IZ=K3 1 ,K32) ,IR= 1 , 1 0) ,DZ ,DR,N,R PROBE 
2 ,V1 ,ZCTR ,FI NT ,FREAL,T 1 ,R K2 
J = 16 
K = 30 
X21:S 
K22:10 
K3 1 = 1 1 
K32=20 
RETURN 
END 


20 
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C INCPSS 

C INCREMENTS TEST POSITION BY STRAIGHT LINE PATH FOR ONE mo 
C 

SUBROUTINE I NCPOS(ZD,X,Y,RD,COSTH,SINTH,VZD,VXD,VYD) 

COMMON/J J/DT /II/RMAX,PI ,IRMAX,VM 

DX=VXD*DT 

DY=VYD*DT 

XrX+DX 

Y=Y+DY 

RD=SQRT(X*X+Y*Y> 

COSTH=X/RD 
SI NTH=Y/RD 
15 ZD=ZD-MED*DT 

RETURN 
END 


C MONTE 

C DISCRETE VARIABLE M/C CHOICE 
C 

FUNCTION MONTE(YID) 
DIMENSION YID (3) 

Y=0. 

DU=0. 

YRtRANDU(DU) 

NY-! 

520 YrY+YID(NY) 

IF(YR.LE.Y)GO TO 600 
NY = NY+ 1 
GO TO 520 
600 MONTE=NY 

RETURN 
END 


C NEVBG7 
C A DTF3 3 

C OFF-LINE DATA PROCESSI NGjLOAD WITH RESNM7 
C 

INTEGER CT,CP,VZ3,VR3,VZ2,VR2 
REAL N 
LOGICAL LSSW 

DIMENSION NI(3P,10)^?2(30,|0>,N3(30,10),M!<30,10>,M2 (30,10) 
2,M3(30,10),SVZ3(23,I'3),VZ3(20,10),SVZ3S(20,10), 

3SVXY3S(20, 10) ,IT3(23,13),YI(3), 

4SVR3 (20, 1 0> , VR3 (20, ! 0) ,M2 0( 1 0) ,M33( 1 0) 

COMMON/I I /R MAX, PI ,IRMAX»VM/JJ/DT /MM/IZMAX»ZMAX/NN/M3 ,SVXY3S, 
2SV73S,SVZ3,SVR3,IZ,IR,N1,N2.N3,M1 ,M2,VZ3,VR3,IT3, 

3 Ml 0(1 0) ,FI NT,FREAL/KX/ZCTR ,RPROBE/LL /SVR2( 1 0, 10), 

4 SVZ2( 10, 1 0) ,VZ2 (10,10) »VR2( 10, 1 0)/III/RK2 ,VMT1 ,UZ2,UR2 

5 ,UZ,UR,V1 , THETA »T3C/JJJ/DR,DZ,T 1 , I ,K,ZPROBE 

6 /KKK/IZC,CT ,CP,S,N/LLL/L1/IJ/FMF?,T1S 
C READ ACCUMULATED PARAMETERS 

REWIND 3 
VRITE(6,50) 

50 FORMAT (15H1 TYPE FMFP,TtS) 

READ (4, 80>FMFP,T 1 S 
80 FORMAT ((El 0.3)) 

VRITE(€,80)FMFP,T1S 


c 

300 WRITE<6,425) 

CALL TAKE (DR,DZ, N> 

VOL =P I* FLOAT (IF:MAX*IRMAX-(IRMAX-1 )* (IRMAX-1 ) ) 

312 DO 340 IZ=1,1ZMAX 

IF(M1(IZ,IRMAX),GT,1) GO TO 33<? 

340 CONTINUE 

Ml<rZMAX f IRMAX)=! 

350 FNORMr N* VOL/FLOAT W Id Z, I RMAX)+I12 (IZ, IRMAX)+M3 Cl c t Lt nflAX)) 

C NEW NUMBER DENSITIES 

DO 401 IR=l,IRMAX 
VOL=PI*FLOATaR*IR-(:H:-l>*(IR-l)) 

DO 360 IZ=!,IZMAX 
RMI=M1(IZ,IR) 

RM2=M2(IZ,IR) 

RM3=M3(IZ,IR) 

N1 (IZ,IR) = FINT*FMGRI' i *RM 1/VOL 
N2(IZ,IR)-FINT*' N0RM*RM2/V0L 
N3 ( I Z , I Fv ) = FI NT* FN OR M*R M3 /VOL 
360 CONTINUE 
C NEW TEMPS & VELS FOR C* ASS 3 
DO 400 lZ=tl,fZMAX 
RM3=M3(IZ,IR)-i 
IZC=IZMAX-IZ+1 
IF(RM3.GT.1.) GO TO 370 
ITZ3=T1 
ITR3=T1 
GO TO 400 

370 VZ3(IZC,IR) = SVZ3(IZC,IR)/RM' 

VR3CIZC,IR) = SVR3(IZC,IR)/RM3 
UZ=VZ3(IZC,IR) 

UR=VR3(IZC,IR) 

ITZ3=ABS((2./RK2)*((SVZ3S(IZC,U)/RM3> 

2-UZ**2) ) 

ITR3=A3S(2.*(SVXY3S(JZC,IR)/RM3«UR*UR)/RX2) 

400 IT3(IZC,IR)= (ITZ3+2*ITR3)/3 
C NEW VELS FOR CLASS 2 

DO 401 IZ=2t,33 
IFCM2(IZ,IR),LE.0) GO TO 401 
IZC=IZMAX-IZ+1 
RM2=M2(IZ,IR) 

VZ2(IZC,I R) =SVZ2 (IZC,IR)/RM2 
VR2CIZC,1R)=SVR2(IZC,IR)/RM2 

401 CONTINUE 
C 

425 F0RMATC18H IF WRITE THIS;SW4) 

R N=FWEAL*FLOAT < N 1 (23,5)+N2(23,5)+N3(23,5))/N 
WRITE (6, 420) RN 
420 FORMAT (2X.F6. 2) 

I F ( . NOT . LSSW ( 4 ) ) GO TO 300 

CALL RESULT (IZMAX.IRMAX) 

GO TO 300 
END 


C NEWBG 

C fCW BACKGROUND DURING SIMULATION 
C 

SUBROUTINE NEWBG(IZMAX,IRMAX) 

INTEGER CT,CP,VZ2,VR2,VZ3,VR3 
REAL N 

COMMON/NtJ/M3(30,l0),SVXY3S(20,10),SVZ3S(20, 10),SVZ3(20,10), 

2 SVR3 <20,10),IZ,IR,N1(30,13),N2<30,10),H3<30,13)» 

3 Ml (30,10),M2C30,10),VZ3(20,10),VR3C20,10),rT3<20,l0) 

4 ,M 10(10), FI NT,FREAL/LL/SVR2 (10,13) 


5 ,SVZ2<10,i0>,VZ2(10 # l&),VR2(10 l I3>/III/RX2,VMTl , 

€ UZ2,UR2 r UZ,UR, VI .THETA ,T3C/I I /RMAX.PI , IRMXD.VM 
7 /KKK/IZC .CT # CP ,S » N/JJJ/DR ,DZ,T 1,1 .X.ZPROBE 
C 

VRITE(6 f 305> 

305 FORMAT (3X.3HTZ3 ,3X,3HTR3) 

VOL=PI*FLOAT(IRMAX*IRMAX-(IRMAX-l>*(IRMAX-n> 

312 DO 340 IZ=i .IZMAX 

IF(Ml(IZ,IRMAX).GT.l) GO TO 350 
340 CONTINUE 

MICIZMAX.IRMAX):! 

330 FNORM= N* VOL /FLOAT (M 1 (IZ »IRMAX)+M2 (IZ f IRMAX)+M3 (IZ.IRMAX)) 

C 

DO 401 IR = 1,IRMAX 

VOL =P I* FLOAT (IR*IR-(IR-l)*(IR-D) 

DO 360 IZ=I, IZMAX 
RM1=M1(IZ,IR> 

RM2=M2(IZ,IR> 

RM3=M3(IZ,IR> 

Nl (IZ,IR)=FINT*FNORM*RM1/VOL 
N2 (I Z , I R > = FI NT*FN0RM*RM2 /VOL 
N3(IZ,IR) = FINT*FN0RM*RM3/V0L 
MI(IZ,IR> = 0 

IF(IZ.L£.20) M2(IZ,IR>=0 
IF(IZ.LE.10> M3(IZ,IR)=0 
360 CONTINUE 

C 

DO 400 IZ = II,IZMAX 
RM3=M3(IZ,IR> 

IZC=IZMAX-IZ+1 

IF(RM3.GT.l.) GO TO 370 

ITZ3=Tl 

ITR3=T1 

GO TO 380 

370 VZ3(IZC,IR>=SVZ3(IZC,IR)/RM3 

VR3(IZC f IR)=SVR3(IZC»IR)/RM3 
UZ=VZ3(IZC,IR) 

UR=VR3(IZC,IR) 

I7Z3=ABS((2./RK2)*((SVZ3S(IZC,IR)/RM3> 

2-UZ**2>)*RM3/(RM3-l.) 

I7R3=ABS(2.*(SVXY3S(IZC,IR>/RM3-UR*UR>/RK2}*RM3/(RM3-I.) 
380 IF(IR.E0.5) VRITEC6,385)ITZ3,ITR3 

385 FORMAT (21 6) 

IT3(IZC,IR)=(ITZ3+2*ITR3>/3 
M3(IZ,IR> = 0 
SVR3(IZC,IR)=0. 

SVZ3(IZC,IR) = 0. 

SVZ3S(IZC,IR)=0. 

400 SVXY3S(IZC,IR>=0. 

C 

DO 401 IZ=21,30 

IF (M2 (IZ t IR).LE.0) GO TO 401 

IZC=IZMAX-IZ+l 

RM2rM2(IZ,IR) 

VZ2(IZC,IR)=SVZ2(IZC,IR)/RM2 
VR2(IZC,IR)=SVR2(IZC,IR)/RM2 
M2(IZ,IR>=0 
SVR2(IZC,IR) = 0. 

SVZ2(IZC,IR) = 0. 

401 CONTINUE 
RETURN 
END 


C NEVCEL6 

C ft XT CELL; TEST FOR SYSTEM BOUNDARY 
C 

SUBROUTINE NEVCEL (Z,X,Y,R ,Vx,VY, VZ.CREF.J) 

INTEGER CT,CP,VZ3,VR3,VZ2,VR2 
REAL N 
LOGICAL LSSW 

DIMENSION N1 (30, ! 0) , N2 (30, 1 0) ,N3 (30, 1 0) ,M 1 (30, 10) ,M2 (30, l 0) 
2,M3(30,13>,SVZ3(20,10),VZ3(20,10),SVZ3S(20,10), 
3SVXY3S(20,1O),IT3(20,10),YI (3), 

ASVR3 (20,10),VR3(20,13) , M2 0(1 3),M30(10) 

COMMON/I I /RMAX,PI, IRMA X,VM/JJ/DT/MM/IZMAX,ZMAX/NN/M3 ,SVXY3S , 
2SVZ3S,SVZ3,SVR3,IZ,IR,N1,N2,N3,M1,M2,VZ3,VR3,IT3, 

3 Mlfl(10),FI NT.FPEAL /XK/ZCTR,R PROBE /LL/SVR2( 10,10), 

A SVZ2(10,I0).VZ2(10.13).VR2(10,10)/III/RK2.VMT1.UZ2,UR2 

5 ,UZ,UR,V1 ,C0STH,T3C,SI NTH/JJJ/DR ,DZ,T l ,1 ,K,ZPKUBt 

6 /KKK/IZC,CT,CP 
C RADIAL BOUNDARY? 

J - 1 

IFCIR.LE.IRMAX) go TO 102 

IF(CT.NE.I) RETURN 

CALL RBOUND (IR ,R ,VX,VY,COSTH,SI NTH) 

C PROBE? 

102 IF(Z.LE.ZPRORE) GO TO 1 08 
IF(R.GT.RPROBE) GO TO 108 
ZRNG=Z-ZCTR 
RADSQ =R *R+ZR NG*ZR NG 

IF(RADSO.GE.(RPROBE*RPROBE)) GO TO 108 
IF (RA NDU(DU) ,GT .CREF) RETURN 
CALL HITSPH(VX,VY,VZ,Z ,R ,X,Y) 

IZ-Z/DZ+l » 

IZC=IZMAX-IZ+1 
IF (CT ,LT #2) CT=2 
Jr3 

RETURN 

C AXIAL BOUNDARY? 

108 IF(IZ.GT.IZMAX) RETURN 

112 IF(IZ.LT.l) RETURN 

IZC=IZMAX-IZ+1 
IF(IZC,GT ,20) GO TO 1010 
C NEXT CELL B/G PARAMETERS 
T3C=IT3(IZC,IR) 

UZ=VZ3(IZC,IR) 

UR=VR3(IZC,IR> 

GO TO 1020 

1010 T3C=IT3(20,IR) 

UZ=VZ3(20,IR) 

UR=VR3(20,IR) 

1020 IF(IZC ,LE. 1 0) GO TO 1120 
UZ2=VZ2(10,IR) 

UR2=VR2(10,1R) 

GO TO 1 125 

1120 UZ2=VZ2(IZC,IR) 

UR2=VR2(IZC,IR) 

1125 J -2 

RETURN 


C NEWPOS 

C POSITION ON ENTRANCE PLANE 

SUBROUTINE NEWPOS (ZD , PD ,THETAD ) 
COMMON/II/RMAX,PI,IRMAX,VM 
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DU=0. 

R2=RANDU(DU) 

ZD=0. 

RD=SQRT(RI)*RMAX 

THETAD=2.*PI*R2 

RETURN 

END 


C NEWV 

C VELOCITY AT ENTRANCE PLANE 

SUBROUTINE NEWVOJZD ,URD ,VX,VY,VZ> 

COMMON/II /RMAX,PI , IRMAX»VM 

DATA A0,Al,A2,A3,A4,A5,A6,A7/-2.55179,4.57444,-2.51282 
2 ,-.92546, 3. 86399, -2. 64662,. 64784, -.02726/ 

DU=0. 

RZ=RANDU(DU> 

C POLYNOMIAL FIT 

X = (AL0G(1./RZ)>**.25 

‘ Y=A0+(Al + (A2+(A3+(A4+(A5+<A6+A7*X>*X>*X)*X)*X>*X>*X 
VT=VM*Y 
VZ=VT+UZD 
C 

RR=RA NDU (DU) 

VRD=URD+VM*SORT (ALOG (I ,/RR) ) 

R THETA :PA NDU (DU) 

THETAD =2.*PI*RTHETA 
VX:VRD*COS(THETAD) 

VY:VRD*SIN(THETAD) 

RETURN 

END 


C PTNR37 

C TEST VELOCITY AFTER COLLISION 

SUBROUTINE PTNR3 (VBZD ,VBRD,COSTHB,SI NTHB,UZD ,URD ,COSTHU ,SI NTHU) 
COMMON/V/VZD ,VXD ,VYD 

C NEW VZ 

VC=0. 

CALL VELS (VC.COSPSI ,THC) 

VSVr 1 . 

CALL VELS (VSW,CSPSIR,THR) 

VSW=0. 

VBX=VBRD*COSTHB 

VBY=V3RD*SINTHB 

VR =SQRT ( (VBZD -VZD >* < VBZD -VZD )+ (VBX-VXD >* (VSX-VXD )+ (V3Y-VYD ) 

2 *<V3Y-VYD )> 

VZD=.5*VZD+.5*VBZD+.5*VR*CSPSIR+VC*C0SPSI 

C NEW VX 

SINPSI=S0RTU.-C0SPSI*C0SPSI) 

SNPSIR=S0RT(1.-CSPSIR*CSPSIR) 

VXD = .5*VXC+.5*V3X+.5*VR*SNPSIR*COS(THR>+VC*SINPSI*COS(THC) 

C NEW VY 

VYD=.5* VYD+.5*V3Y+.5*VR*SNPSIR*SI N(THR)+VC*SI NPSI *SI N(THC) 

RETURN 

END 


C RANDD 

C PSEUDORANDOM NUMBER GENERATOR 


FUNCTION RANDU(DU,X) 

DOUBLE PRECISION X.DXI 
IF(DU.GT.I.) X=DBLE(DU)/1.D5 
X=X*997.D0 
IXI=IDINT(X) 

RXI =FLOA T (I XI ) 

DXI=DBLE(RXI> 

XrX-DXI 

DU=0. 

RANDU=SNGL(X> 

RETURN 

END 


C RBNDS 

C RADIAL SYSTEM BOUNDARY REFLECTION 
C 

SUBROUTI NE RBOUND (IRD,RD,VX,VY,C,S) 

COMMON/II/RMAX,PI,IF!MAX,VM 

IRD=IRMAX 

RD=RMAX 

VDOTRU: VX*C+ VY*S 

VX=VX-2 ,*VDOTRU*C 

VY=VY-2.*VD0TRU*S 

RETURN 

END 


C RESGIV 

C OUTPUT TO DECTAPE 
C 

SUBROUTINE RESULT (IZMA X,I RMAX) 
INTEGER CT,CP 
REAL N 

COMMON/JJJ/DR,DZ,T l/KKK/IZC,CT,CP,S,N 
CALL GIVE(DR,DZ,N) 

RETURN 

END 


C RESWM7 

C OFF-LINE DATA REDUCTION: OUTPUT NORMALIZED PARAMETERS. LOAD 
C WITH NEWDG7 

SUBROUTINE RESULT (11,12) 

INTEGER CT,CP,VZ3,VR3,VZ2,VR2 
REAL N 
LOGICAL LSSW 

DIMENSION N J (30, 1 0) ,N2 (30, 1 0) , N3 (30, 1 0) ,M 1 (30, 1 3) ,M2 <33, 1 0> 
2,M3(30,10),SVZ3(20,1?),VZ3(20, I0),SVZ3S(2O,10), 

3SVXY3S(20, 1 0) ,IT3(20,1 0),YI (3) , 

4SVR3 (20,1O),VR3(2O,!O),M2O(13),M30(10)»RNORM(10),T(10),RN(I3>, 
5 RVZ ( ! 0) t P VR ( 1 0) 

COMMON/I I /RMAX, PI, IPMA X.VM/JJ/DT/MM/IZMAX ,Z*A X/NFJ/M3 ,SVXY3S, 
2SVZ3S ,S VZ3 , S VR3 ,1 Z , I R , N I , N2 , K3 ,M J ,M2 ,VZ3 , VR3 , 1 T3 , 

3 M! 0(10), FI NT,FPEAL/KK/ZCTR,PPRP3E/LL/SV?2(I3,1<?> , 

4 SVZ2( 10,10) ,VZ2 (1 0, 1 0) ,'/R2 ( I 0, 1 0) /III /RED ,VMT 1 ,L'Z2 ,UR2 

5 ,UZ ,UR , VI ,THETA ,T3C/JJJ/DR ,PZ ,T I ,1 ,K,ZPR09E 

€ /KKK/IZC,CT,CP ,S ,M/I J/FMFP ,T IS/LLL/L 1 ,UZM(30> ,URM(30) 
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EQUIVALENCE (KNORMd ) ,R N< 1 > ,P VZ (1 ) ,T ( 1 ) ,R VR ( ! ) ) 

C 

WMEAN<X!,X2,y3,Nl,N2,N3>r(Xl*FL0AT(NI)+X2*FL0AT(N2)+ 

2 X3*FL0AT(N3 ) ) /FLOAT (N 1+N2+N3) 

RK2=5 74. 

DZMFP=DZ/FMFP 
DRKFP=DR/FMFP 
ZCTMFP=ZCTR/FMFP 
WRITE(6, 151) FI NT ,FREAL 
10 FORMAT (1H1,2EI!,2) 

WRITE (6,20) 

23 F0RMAT(18X,8H T TABLE/20X.5HR NORM) 

DO 43 IR=1,I0 

48 RNORM (IR )=DRMFP* FLOAT (IR) 

WRITE (6 ,45 ) (RNORM (IR) ,IR= I , 1 0) 

45 FORMAT (8X,l0F6.i/2X,6H ZNORM) 

C 

DO 100 I Z - 1 1 , 1 ZMAX 
DO 60 IR = 1,I3 
IZC=IZMAX-IZ+I 
RM1=M1 (IZ.IR) 

RM2=M2 (lZ f IR) 

RK3=M3CZ,!R) 

RM=RM1+RM2+RM3 

SVZ1S=RM!*((TI*RK2*(RM1-1.))/(2.*RM1 )+Vl*VI) 

IF(RM.6T.I.) 60 TO 47 

TM=Tt 

GO TO 60 

47 SVRIS=SVZ1S-VI*VI*RM1 

IF«ZC.LE.10> GO TO 50 
UZ2s-VI 
UR2=0. 

GO TO 55 

50 VZ2(IZC,IR) = SVZ2(IZC,IR)/RM2 

VR2<IZC,IR) = SVR2(IZC,IR)/RM2 
UZ2=VZ2CIZC,IR) 

UR2=VR2(IZC,IR) 

55 SVZ2S=RM2* C (T I *RK2*(PM2- 1 .) ) /(2.*RM2 )+UZ2*UZ2) 

SVR2S=SVZ2S-UZ2*UZ2*RM2 
UZM(IZ)= (VI *RM1+UZ2*RM2+SVZ3(IZC,IR) >/RM 

TZ d?M*2.* ( (SVZ I S+SVZ2S+SVZ3SCI ZC , IR) ) /RM-UZMC IZ)*UZM(IZ) )/( (RM- 
2 1 )*FK2 ) 

URM(IZ): OJR2*RM2+SVR3(IZC,IR))/RM 

TR=«M*2.*((SVR1S+SVR2S+SVXY3S(IZC,IR))/RM-URM(IZ)*URM(IZ))/((RM- 
2 I.)*RK2) 

TM=(TZ+2.*TR)/3. 

60 T(IR)=(TM-T1)/(T1S-T1) 

Z NORM=DZMFP*FLOAT (IZ)-ZCTMFP 
WRITE(6,85) ZNORM, (T(IR),IR=I, 10) 
g5 FORMAT(F6.1,2X,|0F6.2) 

100 CONTINUE 

I F (LSSV (0) ) RETURN 
C 

WRITE (6,205) 

205 FORMAT (/) 

WRITE(6,210) 

2! 0 FORMAT ( 1PX,9H RN TABLE) 

DO 270 IZ=!,IZMAX 
DO 220 IR = 1 ,10 

220 RN(IR) = FREAL*FL0AT(N1 (IZ,IR)+N2 (IZ,IR)+N3(IZ,IR) )/N 

Z NORMrDZMFP* FLOAT (IZ)-ZCTMFP 
270 WR I TE ( 6 , 85 ) ZNORM, (RN(IR) ,IR:1 ,10) 

WRITE (6,205 ) 

WRITE (6, 275) 

275 FORMAT (IPX,9H VZ TABLE) 

C 

DO 390 IZ=H,IZMAX 
DO 380 IP=I,I0 
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IZC:IZMAX-IZ + 1 
IFCIZC.LE.10) GO TO 350 
VZ2C=-VI 
GO TO 370 

350 VZ2C=VZ2(IZC,1R) 

370 VZ3C:VZ3(IZC,IR) 

VWEAN=WMEA N(Vl ,VZ2C , VZ3C ,N1(IZ,IR),N2(IZ,IR) 

2 ,N3(IZ,IP)) 

380 RVZ(lP) = 'rEA?J/Vl 

ZRORM=DZMFP* FLOAT (IZ)-ZCTMFP 
390 VRITE(6,85) ZNOPM, (RVZ (IR) ,IR: I , ! 0) 

C 

WRITE (6,395) 

395 FORMAT (//l 8X,9H VR TABLE) 

DO 460 I Z = 1 1 ,1 ZMAX 

DO 450 IR=!,IPMAX 

IZCrIZMAX-IZ+1 

IF(IZC .LE . 10) GO TO 430 

VR2C=0 

GO TO 440 

430 VR2C:VR2(IZC,IR> 

440 VR3C:VP3(IZC,IP) 

VRMN:WMEAN(0.,VR2C,VR3C,Nl(IZ,IR),N2(IZ,IR),N3(IZ,IR)> 
450 RVR(lR):VRMN/Vl 

ZNORM=DZMFP*FLOAT(IZ)-ZCTMFP 
460 WRITE (6,85) ZNORM, (RVR (IR) ,IR: I , 1 0) 

RETURN 


C RESRW7 

C OFF-LINE OUTPUT NON-NORMALIZEO PARAMETERS. 

C LOAD WITH N>:WBG7. 

C 

SUBROUTINE RESULT (IZM,IRM) 

INTEGER VZ3,VR3,VZ2,VR2,CT,CP 
REAL NT , NTOT,N 
DIMENSION NTOT(IO) 

C0MM0N/UN/M3 (30, 1 0) ,SVXY3S(20, 1 0) ,SVZ3S(20, 1 0) ,SVZ3 (20,1 0) 

2 ,SVR3(20,10),IZ,IR,N! (30,10),N2(30,J0),N3(30, 10) 

3 ,M 1 (30, 1 0) ,M2 (30, 1 0) ,VZ3 (23, 1 0) ,VR3(20, ! 0) 

4 ,IT3(20,!3),M10(I0),FINT,FREAL/LL/SVR2(10,I0) 

5 ,SVZ2( 10,t0),V/Z2(I0,10),VR2(10,I0) /JJJ/DR ,DZ,T I /KKK/IZC,CT ,CP 

6 ,S*N/II/RMAX,PI ,IRMAX/MM/IZMAX,Zf1AX/1 1 1 /RK2 ,VMT I,UZ2,UR2 

7 ,UZ,UR, VI , THETA, T3C 
WR1TE(6,402) FI HT.FREAL 

402 FORMAT (6H FI NT: ,E9.2,7H FREAL: ,E9.2) 

IR:5 

WRI TE (6,422) IR 
422 FORMAT (4H IR = ,I3) 

WRITE (6,425) 

425 F0RMAT(3X,3HVZ2,3X,3HVR2,3X,3HVZ3,3X,3HVR3) 

DO 428 IZ=I1,IZMAX 

IZC=IZMAX-IZ+1 

IFdZC.LE.10) GO TO 426 

IVZ2:-Vi 

IVR2:0 

GO TO 428 

426 IVZ2: VZ2 (IZC,IR) 

IVR2=VR2(IZC,IP) 

428 WRITE(6,433) I VZ2 , I VR2.VZ3 (IZC,IR) ,VR3 (I ZC, IR) 

430 FORMAT (416) 

WRITE (6,443) 

440 FORMAT (9H N2 TABLE) 



c wtot 

DO 445 IZ=I,IZMAX 
DO 442 IR=!,IRMAX 
USED HERE TO SAVE STORAGE 

442 

NT0T(IR) = FREAL*FL0AT(H2(IZ,IR)) 

445 

WRITE (6,465 ) (NT OT (IR),IR:I,6) 

455 

WRITE (6, 45 5) 

FORMAT (9H N3 TABLE) 

457 

DO 460 IZrl.IZMAX 
DO 457 IR=I,IRMAX 
8T0T(IR) = FR£AL*FL0AT<N3<IZ,IR)) 

460 

WRITE(6,465)(NT0T(IR),IR=I,6) 

465 

F0RMAT(6E9.2) 

485 

WRITE(6,4B5) 

FORMAT (9H N1 TABLE) 

487 

DO 490 IZ=1,IZMAX 
DO 4R7 I R = 1 ,IRMAX 
NT0T(IR) = FREAL*FL0AT(N1 (IZ,IR)) 

490 

WRITE (6, 465) (NTOT (IR) ,IR=1 ,6) 


RETURN 

END 


C RSTDV7 

C OFF-LINE OUTPUT STANDARD DEVI ATI ON. LOAD WITH STDV7. 

C 

SUBROUTINE RESULT (1 1 ,12 ) 

INTEGER CT,CP,VZ3,VR3,VZ2,VR2 
REAL N 
LOGICAL LSSW 

DIMENSION N1 (3 0, 1 0) ,N2 (30, 1 0) , N3 <30, l 0) ,M 1 (30, 1 0) ,M2 (30, 1 0) 

2, M3 (30, 10) ,SVZ3(20,1 3) ,VZ3 (20, 1 0) ,SVZ3S(20, 10) , 
3SVXY3S(20,!0),rr(20,10) f YI<3), 

4 SVR3 (20, 10), VR 3(20, 10) 

C0MM0N/II/RMAX,PI,IRWAX,VM/JJ/DT/MM/IZMAX,ZMAX/NN/M3,SVXY3S, 
2SVZ3S,SVZ3,SVR3,IZ,IR,N1.N2,N3,M1 ,M2,VZ3,VR3,IT, 

3 Ml 0(1 0) ,FI NT ,FREAL/XK/ZCTR ,RPR03E/LL/SVR2 ( 1 0, < 0) , 

4 SVZ2(10,10),VZ2(10,I0),VR2(I0,10)/III/RK2,VMT1 ,UZ2,UR2 

5 ,UZ,UR ,V1 ,THETA ,T3C/JJJ/DR ,DZ,T 1 ,1 ,K ,ZPROBE 

6 /KKX/IZC,CT,CP,S,N/IJ/FMFP,T1S,XM(30),SSQ(30)/LLL/L1,UZM(30) 

7 URM(30) 

C 

IF(LI-l) 200,50,300 
C 

50 DO 100 IZ=!1,IZMAX 
IZC=IZMAX-IZ+I 
TM=IT(IZC,IR> 

RT= (TM-T 1 )/(T 1 S-T 1 ) 

XM(IZ)=XM(IZ)+RT 
100 SSO(IZ) = SSG(IZ>+RT*RT 

RETURN 
C 

200 DO 270 IZ = 1,IZMAX 

RN=FREAL*FL0AT(N1 (IZ,IR)+N2(IZ,IR)+N3(IZ,IR))/N 
XM(IZ):XM(IZ)+RN 
270 SSO(IZ) = SSQ(IZ)+RN*RN 

RETURN 
C 

300 DO 320 IZ = 11,IZMAX 

IF(L1.GT.2> GO TO 350 
RVZ:UZM(IZ)/V1 
XM(IZ)=XM(IZ)+RVZ 
320 SSQ(IZ)=SSQ(IZ)+RVZ*RVZ 

RETURN 


C 


ouco « «co 
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350 DO 370 1Z = I!,IZMAX 

RVR=URM(IZ)/Vl 
XM(IZ)=XM(IZHRVR 
370 SS0(1Z) = SSQ (IZ)+RVR*RVR 

RETURN 
END 


C SETAK6 

C SET UP INITIAL PARAMETERS. BACKGROUND FROM ANY PREVIOUS 
C ITERATION IF DATA SWITCH 6 SET 
SUBROUTINE SETUP(N) 

INTEGER CT,CP,VZ3,VR3,VZ2,VR2 
REAL N 
LOGICAL LSSW 

DIMENSION N1(30,!0>,N2(30,10),N3(30 # 10),M1(30,10),M2(30,I0) 
2, M3 (30, 1 0) , SVZ3 (20, 1 0) , VZ3 (20, 1 3) ,SVZ3S(20, 1 0) , 
3SVXY3S(20,10>,IT3(20,I0),YI(3) t 
4SVR3 (20,1 0) ,VR3 (2.0, 1 3) ,M2 0( 1 3) ,M3O(10) 

COMMON/I I/RMAX,PI,IRMAX,VM/JJ/DT /MM/IZMAX.ZMAX/NN/MJ ,SVXY3S , 
2SVZ3S,SVZ3,SVR3,IZ f lR,NI,N2,N3,MI,M2,VZ3,VR3,IT3, 

3 MI0(10),FINT,FRFAL/KK/ZCTR,RPROBE/LL/SVH2(10,10J, 

4 SVZ2(I3,1O>,VZ2(10,I3),VR2(!3,10)/III/RK2,VMT1 ,UZ2,UR2 

5 ,UZ,UR ,V1 ,THFTA,T3C/JJJ/DR,DZ,Tl «I ,K,Z PROBE 

6 /KKK/IZC,CT ,CP,S/LLL/L1/KKM/KM 
C 

X:0 

RK2:574. 

C RK2=2*K/M 

S:.I05E-1R 

VMTI=SQRT(RK2*T1> 

FI NT: l.E-17 
FREAL=1.EI7 
PI =3. 14159 

RMAX:DR*FLOAT (IRMAX) 

ZMAX:DZ*FtOAT(IZMAX) 

RPROBE: .005 
ZCTR=ZPROBE+RPROBE 
C 

IF (.NOT. LSSW (6)) GO TO 2 
REWIND 3 
C 

WRITE (6, 100) 

100 FORMAT (13H 4 ITERATIONS) 

READ (4,110) I 
110 FORMAT (13 > 

WRITE(6, 110)1 
DO 120 IT=1,I 

120 CALL TAKE(DR,DZ,N) 

K:KM 

RETURN 


DO 9 IR: 1 , 1 0 
DO S 1Z=1,IZMAX 
IF(Lt.EQ.l) GO TO 3 
Nl(IZ,IR) = f*FINT 
N2(IZ,IR):0 
N3(IZ,Ift':0 
Ml (IZ,IR)=0 
M2(IZ,1R):0 
M3(IZ,IR)=0 
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UU B IL- lm£Yl 

SVZ3(IZ,IR ) = 0 

SVR3(IZ,IR)=0 

SVZ3S(IZ,IR)=0 

SVXY3S(IZ,IR)=0 

IF(LI.EQ.l) GO TO 8 

VZ3<IZ,IR)=3 

VR3(IZ,IR>=0 

IT3(IZ,IR)=TI 

8 CONTINUE 
C 

C 

M10(IR)=0 
DO 9 IZ=!,10 
SVR2(IZ,IR) = 0. ' 
SVZ2(IZ,IR) = 0. 
IF(Ll.EQ.l) GO TO 9 
VZ2(IZ,IR)=-V! 
VR2(IZ,IR> = 0 

9 CONTINUE 
' RETURN 

END 


C STDV7 • 

C A DTF3 3 . ~ 

C OFF-LINE STANDARD DEVlAf I ON & ME AN. LOAD WITH RSTDV7. 

C 

INTEGER CT,CP,VZ3,VR3,VZ2,VR2 
REAL N 
LOGICAL LSSW 

DIMENSION N1 (30, 1 3) , N2 (30, 1 0) , N3 (30, 1 0) ,M1 (30, 1 0) ,M2 (30, ! 3) 

2, M3 (30, 10) ,SVZ3(20.13),VZ3(23,10),SVZ3S(20,I0), 

3SVXY3S(20, 1 0) ,IT3 (20, 1 0) ,YI (3 ) , 

4SVR3 (20, 1 3) , VR3 (20, 1 3) ,M20( 1 0) ,M30 (1 0) 

COMMON/II/RMAX.PI .IRMAX.VM/JJ/DT/MM/IZMAX.ZMAX/NN/M3 .SVXY3S. 
2SVZ3S,SVZ3,SVR3 f IZ,IRT,Nl,N2,N3,MI,M2,VZ3,VR3,IT3 t 

3 Ml 0< 1 0> ,FI HT,FREAL/KX/ZCTR,R PR0BE/LL/SVR2 (10,10), 

4 SVZ2(10,10),VZ2(10, 13),VR2(13,10)/III/RK2,VMT1,UZ2,UR2 

5 ,UZ,UR,V1,THETA,T3C/JJJ/DR,DZ,T1»I,X,ZPR09E 

6 /KKK/IZC.CT ,CP,S,N/LLL/L I ,UZM(30) ,URM(33)/I J/FKFP,T 1S,XM(30), 

7 SSQ (30) 

100 REWIND 3 

VRITE(6,200) 

200 FORMAT ( !4H TYPE L!,IR,NI) 

READ (4,220)L1 ,IRT,NI 
220 FORMAT ((13) ) 

VRITE(S,220)L1,IRT,NI 

VRITE(6,50) ■. o 

50 FORMAT ( I 4H TYPE FMFP*TIS> 

READ (4,80)FMFP,T1S 
80 FORMAT ( (E 10.3)) 

WRITE (6,50)FMFP,T i S 
DO 250 1=1, NI 

250 CALL TAKE(DR,DZ,N) 

ZCTMFP=ZCTR/FMFP 
DZMFP=DZ/FMFP 
DO 290 IZ=!,30 
UZM(IZ)=V1 
URM (IZ)=0. 

XM(IZ)=0. 

290 SSO(IZ)=0. 

C 

DO 410 1=1,10 
300 CALL TAKE(DR,DZ,N> 

VOL=PI* FLOAT (IRMAX*IRMAX-(IRMAX-1 )*(IRMAX-I ) ) 


312 DO 340 IZ=I ,1 ZMAX 

IF CM I (IZ,IRMAX).GT.l ) GO TO 350 
340 CONTINUE 

M1CIZMAX,IRMAX)=! 

350 FNORM=N* VOL/FLOAT CM! CIZ.IRMA X)+M2 CIZ,IRMAX)+M3 CIZ.IRMAX) > 

IR=IRT 

VOL=PI*FLOATCIR*IR*CIR«l)*ClR-I)> 

DO 360 IZ-1 t IZMAX 
RM!=M1CIZ,IR> 

RM2 = M2ClZ i IR) 

RM3=P13CIZ t IR> 

N1 C1Z,IR) = FI NT*FNOR M*RM1 /VOL 
N2 CIZ , IR ) = FI NT*FN0RM*RM2/V0L 
N3CIZ,IR) = FI NT*FN0RM*RM3 /VOL 
360 CONTINUE 

DO 401 I Z ~ 1 1 ,IZMAX 
IZC=IZMAX-IZ+1 
RMIrMl CIZ,IR) 

RM2 = M2CIZ,IR) 

RM3=M3(IZ,IR) 

RMrRMl+RM2+RM3 

SVZIS=RMI*CCT!*RK2*CRMl-i.))/C2.*RMl)+Vl*Vl) 

SVRIS=SVZ1S-VI*VI*RM! 

IFCIZC.LE.10) GO TO 370 

UZ2=-VI 

UR2=0. 

GO TO 375 

370 VZ2 (IZC ,IR)=SVZ2 CIZC, IR) /RM2 

VR2 CIZC, IR)=SVR2 CIZC # IR)/RM2 
UZ2=VZ2CIZC,IR) 

UR2=VR2(IZC IR) 

375 SVZ2S=RM2*C CT i*RX2*(RM2-l .) )/C2«*RM2)+UZ2*UZ2) 

SVR2S=SVZ2S-UZ2*UZ2*RM2 
UZM C I Z ) = C V 1 *RM 1 +UZ 2*RM2+S VZ3 Cl ZC , I R ) ) /RM 

TZrRM*2.*C(SVZ!S+SVZ2S+SVZ3SCIZC,IR))/RM-UZMClZ)*UZMCIZ))/CCRM- 
2 1) *R K2 ) 

URMCIZ)=CUR2*RM2+SVR3CIZC,IR))/RM 

TR=RM*2.*C (SVR1S+SVR2S+SVXY3SCIZC,IR) )/RM»URMCIZ)*URMCIZ) )/C CRM- 
2 I >*RK2) 

IT3CIZC,IR)=CTZ+2.*TR)/3. 

401 CONTINUE 

C 

CALL RESULT CIZMAX.IRMAX) 

410 CONTINUE 

C 

VRITEC6.430) 

430 FORMAT C22H Z NORM MEAN STDV) 

DO 450 IZ = 1,IZMAX 
ZNORM=DZMFP* FLOAT (IZ)-ZCTMFP 
AVG = XMCIZ)/10. 

STDV=SQRTCABSCSSQCIZ)/!0.-AVG*AVG)/9.) 

450 VRITEC6.460) ZNORM,AVG,STDV 

460 F0RMATCF6.1 t 2F6.4) 

GO TO 100 
END 


C TA K£5 

C READ IN ACCUMULATED PARAMETERS FROM DECTAPE. 
C 

SUBROUTINE TAKECDR ,DZ,N) 

INTEGER VR2,VZ2,VR3,VZ3 


REAL N 

COMMON/II/RMAX.PI ,IRMAX,VM/MM/IZMAX,ZMAX/NN/M3 C30,l 0) , 
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2 SVXY3S(20, 1 0) ,SVZ3S (20, 1 0) ,SVZ3 (20, 1 0) ,SVR3 (20, 1 0) , 

3 IZ,IRT,N1 (30,10),N2(30,l0),N3(30,10),Ml (30, 1 0) ,1*12 (30 ,1 0) , 

4 VZ3(20,!0),VR3(20,10),ITN3(20,10), 

5 M10(10),FINT,FREAL/LL/SVR2(10,l0),SVZ2(!iM0), 

6 VZ2(13, 10) ,VP2(10,I0)/KK/ZCTR,RPROBE/III/RK2, 

7 VMTI ,UZ2,UR2 ,UZ ,UR ,V1 , THETA, T3C 

8 /JJJ/DRD,,DZD ,T 1 
VRITE(6, 10) 

10 FORMAT (10H DT3;SW 17) 

11 CONTINUE 

IF(. NOT .LSSWd 7) )G0 TO 1 1 

J = 1 

X = 13 

K21 = l 

K22=5 

K3 1-1 

K32=10 

0020 L= 1 ,2 

READ (3)((M1(IZ,IR),IZ = J*K),IR=1»10) ,IZMAX 
READ(3)((M2(IZ,IR),IZ=J,X),IR=I,10) 

READ (3 ) ( (M3 ( IZ, I R) ,1 Z= J ,K) ,1 R: 1 , 1 0) ,1 RMA X 

READ (3) ( (SVR2 (IZ,IR),SVZ2(IZ,IR),IZ=K21 ,K22) ,IR= I , 1 0) ,PI 

READ(3)((SVR3(IZ,IR),IZ=K31,K32),IR=I ,10),RMAX 

READ (3) ( (SVZ3 (IZ,IR),IZ;K3I ,K32) ,IR= 1 ,10) 

READ (3)( (SVZ3S(IZ,IR) ,IZ = X3 1 ,K32),IR:1,10),ZMAX 
READ (3) ( (SVXY3S(IZ,IR) ,IZ=K3I,K32),IR-1,10),OZ,DR,N,RPROBE 
2 ,V1 ,ZCTR ,FI NT ,FREAL ,T 1 ,RK2 
J = 16 
K=30 
K21 = S 
K22= 1 0 
K3 1 = 11 
20 K32=?0 

RETURN 
END 


C VELS7 

C RANDOM VELOCITY VECTOR WITH SPHERICALLY SYMMETRIC MAXWELLIAN 
C DISTRIBUTION. 

SUBROUTINE VELS(VCD,COSPSI,THCD ) 

C OMMON/II /RMA X, PI ,IRMAX,VM 
C VCD : l {UNIT VECTOR ONLY. 

IF(VCD.EQ.I.) GO TO 10 
R7=RANDU(DU) 

R8=RANDU (DU) 

ALG=AL0G(1./(R7*RB)) 

VCD=VM*SQRT(ALG) 

10 R9=RANDU(DU) 

COSPSI=I ,-2.*R9 
R10=RANDU(DU) 

. THCD =2.*PI*R 1 0 
RETURN 
END 



C VREL7 

C MEAN EFFECTIVE COLLISION SPEED. 

C 

FUNCTION VREL(VZA,VXA,VYA,VZB,VRB,COSTHB,SINTHB) 

DIMENSION PSI (20) 

COMMON/II/RMAX,PI,IRMAX,VM 

DATA PSI (1), PSI (2), PSI (3) ,PSI ( 4) ,PSI (5) ,PSI (6), 

2 PSI (7> ,PSI (8),PSI (9) ,PSI ( 1 0) .PSI (11) ,PSI ( 12) , 

3 PSI (13), PSI (14), PSI (15)/. 2007,. 4053,. 61 78,. 8420, 1.0813, 

4 1.3390,1 .6182,1.9213,2.2507,2.6083,2.9958,3.4145, 

5 3.8654,4.3494,4.8671/ 

C 

X6=4.*VM*VM/1 .771 
CX=VXA-VRB*COSTHB 
CY=VYA-VRR*SI NTHB 
CZ=VZA-VZB 

C=SQRT(CX*CX+CY*CY+CZ*CZ) 

S=C/VM 

C 

IF (S.LT .1.5) GO TO 5 
ERFzl. 

GO TO 17 

5 IF(S.GE .0.1) GO TO 7 

ERF=S 
GO TO 17 
C 

7 KK=INT((S+.05)*10.) 

PS:PSI(KK) 

GO TO 20 
C 

17 PS=S*EXP(-S*S)+(2.*S*S+l.)*.885*ERF 

20 VREL=PS*X6/C 

RETURN 
END 


